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1 Introduction 


These notes should be used to support the lectures for B5.2: Applied Partial Differential 
Equations. First, however, we pause to review some basic concepts relating to partial differ- 
ential equations and to remember some basic methods for solving ODEs and PDEs. Some 
worked examples are also included to refresh your memory. 


What is a PDE? 


A partial differential equation (PDE) is an equation for some quantity u (the dependent 
variable) which depends on two (or more) independent variables and involves partial deriva- 
tives of u with respect to at least some of these independent variables. In this course, we 
limit attention to just two independent variables, so u = u(x, y), and at most two derivatives. 
Then our PDE can be written in general terms as follows: 


ED Uy hg Uys, Wares Wingy: Uy ) = 0. 


Notes 


1. In applications, x, y often represent spatial variables and a solution may be required in 
some region 22. In this case, there will be some conditions to be satisfied on the domain 
boundary OQ; these are called boundary conditions (BCs). 


2. In applications, one of the independent variables may represent time, t say. Then there 
will be initial conditions (ICs) to be satisfied (i.e., u may be specified everywhere in 
Oak b=). 
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3. In applications, systems of partial differential equations (PDEs) can arise involv- 
ing the dependent variables uy, u2,...,Um, m > 1 with some (at least) of the PDEs 
involving more than one dependent variable. 


4. It is common to use subscripts as shorthand for partial derivatives, 1.e. 
Ou Oru 
ao jie = 1 
Ox’ wy OyOx’ (1) 
5. We will assume that wu is sufficiently smooth for all the required partial derivatives to 
exist and to be independent of the order in which the derivatives are performed, i.e. 
Oru Oru 
— = — Uny = Uyx: 2 
OyOx — OxOy’ ve = (2) 
Definitions 


e The order of the PDE is the order of the highest (partial) differential coefficient in the 
equation. 


As with ODEs, it is important to be able to distinguish between linear and nonlinear equa- 
tions: 


e A linear equation is one in which the equation and any BCs or ICs do not include any 
product of the dependent variables or their derivatives; an equation that is not linear is 


nonlinear. 
Ou Ou : ; : 
Bt + Cae 0 first order linear PDE (the simplest wave equation), 
ee 
ui Ou é : . 
— += 3 = ®(zx,y) second order linear PDE (Poisson’s equation). 
Ox? Oy? 


e A nonlinear equation is semilinear if the coefficients of the highest derivative are 
functions of the independent variables only. 


(a+ jot Ou _ 3 
w+ a)a ee 
Pu, ce iam Ou ,0u_ 4 
Toq2 | UT Wa On ee 


e A nonlinear PDE is quasilinear if it is linear in the derivatives of order m, with 
coefficients depending only on its independent variables, x, y say, and partial derivatives 
of order < m. 


du\*| du Oudu O?u du\*| du 
1+ 2 + }14 = 0. 
Oy Ox? Ox Oy OxOy Ox Oy? 
e Principle of superposition. A linear equation has the useful property that if uv, and 
uz both satisfy the equation then so does au; + Bug for any a, € R. This is often 


used to construct solutions to linear equations (for example to satisfy BCs and ICs, cf 
Fourier series methods). This is NOT true for nonlinear equations. 
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Examples 


Wave Equations Waves on a string, sound waves, waves on stretched membranes, electro- 
magnetic waves, etc 


Oru ait Ou 
Ox2 2. OR?” 


or, more generally, 


where c is a constant, the wave speed. 


Diffusion or Heat Conduction Equations 


au _ Ou 
Ot Ox?’ 
or, more generally, 
Ou 9 
a> KV~U, 
or even 
Ou 
a V -(KVu) 


where & is typically a constant (diffusion coefficient or thermal conductivity). 


Both the wave equation and the diffusion equation are linear equations that involve 
time ¢ as an independent variable. They require initial and boundary conditions for 
their solution. 


A physical interpretation of a function of two variables might, for example, relate to 
the temperature in a metal bar and how this changes over time. Imagine that one end 
of a a metal bar is placed next to a heat source. In order to describe the response to 
heating the bar, we need to know: (a) the temperature at a given position, and (b) the 
temperature at a given time. In other words, the temperature T in the bar should be 
viewed as a function of two real variables, position x and time t, say. From this, we 
can draw a better understanding of the meaning of a derivative of a function of two 
variables. For example, we can describe 7; as ”the rate of change of temperature over 
time at a fixed position of x”. 


Laplace’s Equation Another example of a second order linear equation is Laplace’s equa- 


tion: 2 92 
u u 
—++ = 5 = 0. 3 
Ox? . Oy? (3) 
or more generally 
V2u = 0. 


This equation usually describes steady processes and is solved subject to BCs. 
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Other Common Second Order Linear PDEs 


e Poisson’s equation is simply Laplace’s equation (homogeneous) with a known source 
term (e.g. electric potential in the presence of a charge density): 


V-u= ©. 


e The Helmholtz equation may be regarded as a stationary wave equation: 


V2u t+ k?u =0. 


e The Schrodinger equation is the fundamental equation of physics for describing quan- 
tum mechanical dynamics. The Schr6dinger wave equation is a PDE that describes how 
the wavefunction of a physical system evolves over time: 


Ou 
-V7ut+tVu=i-. 

UutVu=2t DE 
Nonlinear PDEs 


e An example of a nonlinear equation is Fisher’s equation for the propagation of reaction- 
diffusion waves: 


0 oO? 
FT = 5 tu(l—u) 2nd order PDE. 
e The following PDE describes nonlinear wave propagation 


Ou Ou 
Ot + (ut c) a. =0 Ist order PDE. 


This is an example of a Conservation equation 
ut F(u)e = 8, 
where F'(u) is known as a flux function for some density u(x,t) and s is a source term. 


e An example of a quasilinear equation is 


e The following PDE is semilinear 


Ou | 3, 
U5, + 4 
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Systems of PDEs 


e Maxwell equations constitute a system of linear PDEs: 


p ., Loe 

V-E=* Vx B= aca 

ea x HIT 3p 
V-B=0, yee 
Ot 


In empty space (free of charges and currents), this system can be rearranged to give the 
equations of propagation of the electromagnetic field 


2B 
Ot? 


OE 


pe = CeV7B. 


= VE, 

e Incompressible magnetohydrodynamics (MHD) equations combine Navier-Stokes 
equation (including the Lorenz force), the induction equation and the solenoidal con- 
straints: 


“7 +U-VU =—-VII+B-VB+eU+F, 
OB 
a = VX (Ux B)+nV°B, 


V-U=0, V-B=0. 


e Both of the above systems involve partial derivatives in space and time and, therefore, 
require ICs and BCs for their solution 


When can a PDE be solved? 


Having defined a PDE, the general idea is to solve for the function u(x, y). However, a PDE 
does not in general determine u uniquely by itself — boundary conditions are also required. 
A PDE supplemented by appropriate boundary conditions is often referred to as a system or 
a problem. For example, a simple problem involving Laplace’s equation is 


Ou Ou 
oe oR y>O0 


u— 0 y —> Oo. 


We would like to determine a unique solution to the above problem. We have the PDE 
and boundary conditions, but is this sufficient information to achieve our aim? Maybe these 
particular boundary condition will lead to no solutions or multiple solutions, or maybe the 
solutions will be extremely sensitive to small changes in the boundary conditions. It is of 
pointless to attempt to solve a problem for which no solution exists, that many solutions 
exist or to find unusual sensitivity of the solution to the boundary conditions. 

Faced with a problem like (4), the most important question to ask is: is it well posed? 

To be well posed, a problem must have the following three properties: 


1. a solution u(x, y) exists; 
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2. the solution is unique; 
3. the solution depends continuously on the boundary data. 


The first of these is obvious: there is no point in trying to find a solution that does not exist. 
If a problem is physically motivated, and wu represents a physical quantity, then we would 
expect u to have a unique well-defined value at each point. If it does not, it suggests that a 
boundary condition or other constraint is missing from the problem. 

To illustrate the final condition, suppose we vary the function f(x) in (4) by a small 
amount and ask whether the corresponding variation in the solution is similarly small. If it 
is not, then the solution of the problem is practically impossible, since any numerical errors 
in f(x), however small, lead to large errors in the solution. 


Suggested textbooks 


J Ockendon, S Howison, A Lacey & A Movchan, 1999 Applied Partial Differential 
Equations. Oxford. 


R. Courant & D. Hilbert, 1989 Methods of Mathematical Physics Volume II. Wiley. 
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2 First order quasilinear equations 


2.1 Definitions 


In this section we consider partial differential equations (PDEs) of the following form 
O 0 
a(x,ysu) a + B(x, ¥su)a = (2,95). (5) 


Here x and y are independent variables, a, 6 and c are given smooth (i.e. continuously 
differentiable) functions, and u(z,y) is a scalar function for which we would like to solve. 
Equation (5) is known as a first-order quasilinear partial differential equation: first-order 
since there are no second or higher derivatives and quasilinear because it is linear in its 
highest derivatives (there are no nonlinear terms like (Qu/0xr)”). Equations of this type arise 
in many areas of mathematical modelling, including fluid mechanics and traffic flow. They 
also provide a relatively straightforward introduction to some important concepts, such as 
Cauchy data, characteristics and weak solutions, that will be applied to more complicated 
equations later in the course. 

Two special cases of (5) are worth mentioning. First, if a and b are independent of u, then 
(5) becomes the semilinear equation 


Ou Ou 
atst) 5 OU) a = c(a, y, 2). (6) 


If it also happens that c is a linear function of u, then we have 


anu) 5e + On) Se = a(e.yu + Gle.9), (7 
which is a linear equation. It is generally the case that linear equations are significantly 
better-behaved and easier to solve than nonlinear ones. 

Equations like (5) often describe the evolution of a quantity wu (representing e.g. traffic 
density or fluid velocity) in space and time. In such cases, to emphasise the fact that one 
independent variable represents time, we will use « and t as independent variables instead of 
x and y, writing (5) as 


a(t, 1) + bla, tu) = ca, 6,0): (8) 


2.2 Characteristics 
2.3. Geometric definition 


We can think of the solution u(z,y) we are seeking as defining a surface z = u(z,y) in 
three-dimensional space. The normal to this surface is in the direction 


nx V(u(z,y) — 2) = du (9) 


B5.2 Applied Partial Differential Equations 9 


Characteristics 


(a,b,c)" 


Initial . 
curve 
a7, Ee — x 
ae Characteristic 
projections 


y 


Figure 1: Schematic showing the characteristics, parameterised by 7 and pointing in the di- 
rection (a, b,c)", emerging from the initial curve, which is parameterised by s. The projection 
of the initial curve onto the (x,y) plane is P and the projection of the characteristics onto the 
(x,y) plane are the characteristic projections. 


and the PDE (5) can therefore be written as 


b i= 0, (10) 


It follows that the vector (a,b,c)’ is everywhere tangent to the solution surface. 
We can construct curves that are everywhere tangent to (a,b,c)’ by solving the simulta- 
neous ODEs 
dz dy du 
—=alzr,y,uU), —— = d(z,y, u), 7 
dr (2, y,u) dr (2) 95%) dr 
Such curves are called characteristics of the PDE (5). Their projections onto the (x, y) plane, 
i.e. the plane curves (x(7),y(7)) are called characteristic projections. 


= CF; yy it). (11) 


Solution by characteristics 


Suppose, as before, that u is specified along some curve I in the (x, y) plane, i.e. we are given 
u = uo(s) when x = xo(s) and y = yo(s) where s parameterises [. As shown in Figure 1, 
this defines a initial curve in three-dimensional space, through which we require our solution 
surface to pass. For any fixed value of s, we can find a characteristic that passes through the 
point (x0(s), yo(s);, wo(s))* (and is everywhere tangent to (a,b,c)") by solving (11) with the 
initial conditions 


© = 29(5); y = yo(s), U = U0(s) ab =O, (12) 


As s is varied, these characteristics sweep out the desired solution surface. Put another way, 
the initial-value problem (11, 12) determines in principle three scalar functions x(s,7), y(s, 7) 
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and u(s,7), and the vector 
) 
r= | y(s,7T) (13) 
) 


defines the solution surface, parametrised by s and rT. 

The theory that exists for ODEs may be applied directly to the system (11). For example, 
Picard’s theorem tells us that there is a unique solution to (11) satisfying the initial conditions 
(12) provided the right-hand side (a,b,c)” is bounded and satisfies a Lipschitz condition in 
u. However, we can easily construct examples for which these conditions fail and the solution 
either blows up or becomes nonunique at some distance from the initial curve. 

Equation (13) is called the parametric form of the solution. It may be possible to eliminate 
s and 7 from (13) to obtain the solution surface in the implicit form G(x, y,u) = 0. Finally, 
if this implicit equation can be solved for u, then we obtain the solution in the explicit form 
u = u(z,y). Explicit solutions are the most convenient, but are often impossible to obtain in 
terms of elementary functions. 


Example 1 Consider the PDE 


Ou Ou 
oe ee 14 
rae ee (14) 
subject to the boundary data u=0 whenx+y=0. The characteristics satisfy 
dx dy du 
dr dr dr , ) 
and the boundary data lead to the initial conditions 
L=8, y=—-S, u=0 at T = 0. (16) 
Hence we find the parametric solution 
L=st+T, y=—s+rT, U=T, (17) 
and it is straightforward in this case to eliminate s and 7 to obtain the explicit solution 
u+y 
= —_*, 18 
: (18) 


In Example 1, the PDE (14) is semilinear. In such cases, the characteristic projections 
satisfy the ODEs 


dx dy _ 
dr = a(x, y), dr = b(x, y), (19) 


which are independent of the solution u. The standard theory of phase planes may be applied 
to the ODEs (19); for example, there is in general a unique characteristic projection through 
each point in the (x,y) plane except at critical points where a and b are both zero. Once 
(19) have been solved to find the characteristic projections in the (x,y) plane, we find that u 
satisfies the decoupled ODE 
du _ 
a 
along each characteristic projection. 
For general quasilinear equations, the characteristic projections depend on the solution; 
the three ODEs (11) are coupled and must be solved simultaneously. 


c(x(r), UT). u) (20) 
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Example 2 Solve the PDE 


Ou Ou 
uae =] 21 
a ae oe 
for u(a,t) int > 0, subject to the initial condition u =a ont =0. 
The characteristics are given by 
dt daz du 
drt , dr m dr ; ( ) 
and the initial data may be parametrised by 
t=0, r= Ss, u=s at T = 0. (23) 


Solving for t first, we see that t = 7 and thus we may replace T by t henceforth. The initial-value 
problem for u has the solution 


u=st+t, (24) 
so that the problem for x becomes 
d 
W=s+t, x=s whent=0, (25) 
dt 
whose solution is 
r=s+st+3t. (26) 


Now we can solve (26) for s and substitute it into (24) to obtain the solution in explicit form: 


a+t+ 3? 
aoe, 27 
OT +t oie 


Alternative method of solution 
The characteristic equations (11) may be rearranged to give 


dx dy du 
— — » 2 
Gia te au oo 


Suppose we can spot two linearly independent first integrals of these ODEs, of the form 
f(x,y,u) = const and g(x,y,u) = const. Then the general solution of the PDE (5) may be 
written in the implicit form 


f(x,y, u) = F(g(#,y,u)), (29) 


where F is an arbitrary function. 


Example 3 Return to the problem considered in Example 2. The characteristic equations may be 
written as 


dt dx _ du 


= 30 
1 U 1 oY 
and then rearranged to two ODEs: 
du dx 
dt” di ey 


These may be integrated to give 


u=t+C,, c= 5t?+Cit+Cr, (32) 
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where C, and C2 are constants. Our two first integrals are, therefore, C) = f(a,t,u) = u—t and 
C2 = g(2,t,u) = x—3t? Cit =x— 40? —(u—t)t. The general solution is found by setting f = F(g), 
which leads to 
u=t+ F (a+ $i? — ut), (33) 

where F is an arbitrary function. It may readily be verified that any u(x,t) satisfying the implicit 
equation (33) is a solution of (21). 

The function F is found by fitting the initial data: u = x when t = 0 leads to F(x) = a, that is 
u=t+ot+ st? —ut, which reproduces the solution (27). 


This procedure works because the equation f(x,y,u) = const defines a one-parameter 
family of surfaces, as does g(x,y, u) = const, and characteristics are lines of intersection be- 
tween one member from each of these two families. Now, any surface defined by an equation 
of the form f = F(g) has the property that f is constant whenever g is constant. It follows 
that such a surface is composed of a family of characteristics, as indicated in Figure 1, and is 
thus a solution surface for the PDE (5). 


Example 4 For the PDE 


yus — ru a L-Y, (34) 
the characteristic equations 
da dy du 
as = yu, oe = —ZUu, ae =z2—-Y; (35) 
may be rearranged to give [Exercise] 
= (a? +y") = . (u? + 2x + 2y) =0. (36) 
dr dr 
It follows that the general solution is 
u? = —22 — 2y + F(z? +7), (37) 


where F is an arbitrary function. 


2.4 Cauchy data 
Geometric interpretation 


The term Cauchy data refers to the boundary data that, when applied to a PDE, in principle 
determine the solution, at least locally. For the first-order quasilinear PDE (5), Cauchy data 
is the prescription of u on some curve I in the (x,y) plane, that is we set wu = uo(s) when 
x = xo(s) and y = yo(s) where s parametrises [. The combination of the PDE (5) and 
Cauchy data is called the Cauchy problem. For the moment we assume that xg, yo and uo are 
smooth functions of s (although there are interesting cases where this is not true, e.g. where 
I’ has corners) and that there are no values of s for which x25(s) = yo(s) = 0 (this ensures 
that s is a sensible parameter for I). 

We have seen that the method of characteristics, outlined in section 2.2, usually allows a 
solution surface to be constructed in a neighbourhood of I’. However, the procedure fails if 
it happens that the initial curve I is at any point tangent to (a,b)’. If this occurs, then the 
characteristic projection, instead of propagating away from the initial curve, points along it. 
Thus the data uo(s) given on I’ will not in general agree with the ODE (20) satisfied by u in 
the direction (a, b). 
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Example 5 We return to the PDE (14) from Example 1, namely 


Ou Ou 
sales os ie | 38 
Ox a Oy , co 
whose general solution is [Exercise] 
u= 24+ Fwy) (39) 


where F is an arbitrary function. 
Now we attempt to fit three different sets of initial data and thus determine the function F. 


1. u=0 onz+y=0 
This is the case considered in Example 1, in which the initial curve is normal to the characteristic 
projections. The initial data gives F =0, so the solution (18) is reproduced. 


2 u=O0onzrz=y 
This time the initial curve is a characteristic projection. When we attempt to fit the initial data 
we find F(0) = —2, which is impossible. In this case there is no solution. 

3 U=Xonnx=y 
Again, the initial curve is a characteristic projection, but this time we have x = x + F(0) so 
that F can be virtually anything so long as F(0) = 0. This is the nongeneric case in which it 
just happens that the initial data and the characteristic equation (20) agree, and the solution is 
consequently nonunique. 


Example 5 illustrates the following three possibilities. 


1. If is not tangent to a characteristic projection, then there should be a unique solution, 
at least locally. 


2. If [is at any point tangent to a characteristic projection, then there is in general no 
solution. 


3. There is, however, an exceptional case in which the data for u specified on T agree with 
the ODE (20) satisfied by u along characteristic projections. If this happens then there 
is a nonunique solution. 


Cauchy—Kowalevski theorem 


A necessary condition for a unique solution u to exist in a neighbourhood of I is for the first 
derivatives of u to be determined on [’. Differentiation of ug and use of the chain rule leads 


to 
dug _Oudzg | Oudyo 


ds Oxds dy ds’ 
The partial differential equation (5) and (40) form a pair of simultaneous equations for 0u/Ox 
and Ou/Oy on the curve I. We can therefore solve uniquely for these first derivatives so long 
as the determinant of the system is nonzero, i.e. 


(40) 


a b a a 
day dyo| =a i. b 7 #0. (41) 
ds ds 


If this condition is satisfied, then both u and its first derivatives are uniquely determined on 
the curve I’, which is clearly the first step in extending the solution away from I’. Notice that 
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the criterion (41) is equivalent to requiring I not to be tangent to a characteristic projection, 
as argued above via geometrical reasoning. 

When the determinant in (41) is zero, there is either no solution for Ou/Ox and Ou/Oy 
or an infinite number of solutions (this is an instance of the Fredholm Alternative). By 
eliminating between (5) and (40) we find that 


Ldap _ ldyo _, 1duo . 
7 l 42 
a ds b ds * c ds ae no solution, (42a) 


1 dxo 1 dyo 7" 1 duo 


ads bds cds 


=> many solutions. (42b) 


The latter equality is the exceptional case, seen in Example 1, in which the variation of u along 
I’ just happens to agree with the differential equation (20) satisfied along the characteristic 
projection. 

The process outlined above can be continued to obtain higher derivatives of u. If Ou/Ox 
is known, for example, then further differentiation with respect to s gives 


d /0u O07u dx O?u dyo 
= 4 
ds (5) Oe ds ” OxOy ds’ =) 
while differentiation of (5) with respect to x yields 
aot Oru _Oadu , Obdu | Oa (du 2 Ab Ou Ou _ Oc, Oc Ou (44) 
Ox? OxOy | Oxdx Oxdy Ju \ Ax dudx Oy Ox Oudx 


Now we have a pair of simultaneous equations for 02u/Ox? and 0?u/Oxdy. The condition for 
this system to have a unique solution is identical to (41). 

So long as a, 6 and c are analytic, so that this differentiation may be continued indefinitely, 
we can continue this argument to show that the condition (41) allows the derivatives of u to 
all orders to be defined uniquely at [. Thus a Taylor series for u(x, y) may be constructed 
about the initial data curve T, and it can be shown that this series has a nonzero radius 
of convergence. This is the starting point for the proof of the Cauchy—Kowalevski theorem, 
which states that (5) has a unique analytic solution in some neighbourhood of I’, provided a, 
b and c are analytic and satisfy the condition (41). 


2.5 Domain of definition 
Bounded initial curve 


In Example 2, we are given u = x along the whole z-axis. In general, however, the initial 
data may only be given on a finite or semi-infinite initial curve T’. In such cases, the solution 
is only defined in the region penetrated by characteristic projections that intersect I. This 
region, which is bounded by the characteristic projections that pass through the end points 
of I’, is called the domain of definition. 


Example 6 Solve the partial differential equation 


Ou Ou 
DE + tun =u (45) 


for u(x,t) int > 0, subject to the initial condition u= a whent=0,0<a<1. 
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Figure 2: The characteristic projections given by equation (47). 


The characteristics are given by 


dt dx du 
—_— = 1 — = —— 
dr ; Go a 29) 


and the initial data may be parameterised by t 0, x 5, U s,0<s <1. The solution in 
parametric form is 


x = sexp(s(e’ — 1)), u= se’, O<s<1, (47) 


and the characteristic projections are shown in Figure 2. The domain of definition is the region 
0 <a <exp(e! —1). 
Note that s may be eliminated from (47) to obtain the solution in the implicit form 


a = uexp(u—t—ue ‘), (48) 


but there is no explicit formula for u(x,t) in terms of elementary functions. 


Blow-up 


The domain in which the solution is defined may be further restricted if the solution develops 
a singularity, such that the PDE (5) ceases to make sense. For example, u may blow up a 
finite distance from the initial curve [. The method of characteristics reduces the partial 
differential equation (5) to the system (11) of ordinary differential equations. Since nonlinear 
ODEs may certainly give rise to solutions that blow up, the same is true of nonlinear PDEs, 
even those that are semilinear. 


Example 7 Consider the equation 
Ou Ou 3 


ey ee 49 
oat oe (49) 
subject tou=y onx=0,0<y<3. The characteristic equations 
d d 
= ei of =] ow (50) 


dr” dr . dr : 
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Domain of 
1| definition ; 


OF Bia edi giahe -hanavah seas eal reed al al Ge bua. fer aca Apes Rial Sie Gaia Boa Rene, 

-1 sl . [ 1 1 

-1.0 -0.5 0.0 0.5 1.0 
wb 


Figure 3: Domain of definition for Example 7. 


and initial datax =0, y= 8s, u=s onT=0,0<85< 8, have the solution 
8 
L=T, =s4T, uu = —_., 0<s<3. 51 
¢ V1 — 287 oo 
We may solve explicitly for u to obtain 
= cas (52) 
1 — 2a(y — x)? 
This blows up on the line s = 1/V2r, i.e. on the line 
+ : (53) 
=x¢+———. 
¥ V 2x 


The domain of definition, bounded by this curve and the characteristic projections y = x andy = £+3, 
is ulustrated in Figure 3. 


Nonuniqueness 


The domain in which the solution is properly defined may also be limited by u ceasing to be 
a unique function of x and y. Provided the coefficients a, b and c are well-behaved and u 
does not blow up, the method of characteristics outlined in section 2.2 always allows us to 
determine in principle the solution in parametric form: (x(s, T), y(s,7), u(s,7)). Then u 
may in principle be found as a function of x and y so long as there is a unique transformation 
from (s,7) to (x,y). By the Inverse Function Theorem, a sufficient condition is that the 
Jacobian of the transformation be finite and nonzero: 


Ox Ox 
Or 0s Oy Ox 
— — ae aly ee ——, 4 
es By dy an, ba, #0, ore) (54) 
Or Os 


Note that this reproduces the condition (41) for u to be determined in the neighbourhood of 


r. 
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Figure 4: The characteristic projections for Example 8. 


A unique correspondence between (s,7) and (2,y) implies that a unique characteristic 
projection passes through each point in the (x, y) plane. Where J becomes zero, this typically 
signals that the characteristic projections start to cross each other. For semilinear equations, 
this can only happen at critical points of the phase-plane problem (19). 


Example 8 Solve the PDE problem 


Ou Ou 


r— +y— =), u=yonx=1,0<y<1, (55) 
Ox Oy 


and determine the region in which the solution is defined by the boundary data. 
The characteristic equations 


dx dy du 
ao = ———— 56 
dr dr” dr 0, oo 
and initial datax =1, y=s, u=s, att =0,0< 85 <1, lead to 
r=e’, y = se’, U=s, O0<s<l. (57) 


We can eliminate s and 7 to obtain the explicit solution u = y/x in0 < y/a <1. This solution is 
evidently not uniquely defined at the origin where, as shown in Figure 4, the characteristic projections 
all cross and where J becomes zero. We cannot continue the solution beyond this point, so the domain 
of definition is0<y/a<1,4>0. 


For more general quasilinear equations, the characteristic projections depend on the so- 
lution u, so the restriction that they may only cross at critical points no longer holds. The 
generic situation is that J = 0 along curves in the (x,y) plane. On these curves, the solution 
surface starts to fold over itself such that u ceases to be a single-valued function of x and 
y. Since u usually represents a physical quantity (such as pressure, temperature or asset 
price), it cannot be multivalued. Moreover when the solution surface develops a fold, the first 
derivatives of u become unbounded, so the PDE (5) ceases to make sense. For these reasons, 
we have to cut off the domain of definition along any curves on which J is zero. 
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Figure 5: (a) The characteristic projections for Example 9. (b) The domain of definition, 
bounded by the curve on which J = 0. 


Example 9 Solve the PDE problem 


Ou Ou 

— —=0 t>0 58 

ap ag ; (58) 
u =sin(2), O<a<27, t=0, (59) 


and find the region in the (x,t) plane where the solution is uniquely defined. 

The solution is u = sins, t = T, = s+7Tsins in parametric form, or u = sin(x — tu) in 
implicit form. The characteristic projections (found by fixing s and varying T) are straight lines and 
are illustrated in Figure 5(a). We can see that they start to cross a finite distance from the initial data 
t=0. The Jacobian is O(a, y)/O(s,7T) = 1+tcoss. The curve on which J = 0 is, therefore, given 
parametrically by x = s—tan(s), t = —sec(s) and is illustrated in Figure 5(b). The solution is defined 
in the region bounded by this curve, the characteristic projections x = 0 and x = 27, andt =0. 

In Figure 6, we visualise the solution by plotting snapshots of u versus x at different times t. The 
initial u = sin(x) steepens as t increases from zero, becoming multi-valued for t > 1. When t = 1, the 
Jacobian first reaches zero at x = 1, where |Ou/Ox| becomes unbounded. 


A curve on which J = 0 may also be viewed as an envelope of the characteristic projections. 
Given a family of curves F'(z,y;A) = 0, where \ is a scalar parameter, an envelope is a 
curve that at each point meets one of the family tangentially, and is determined from the 
simultaneous equations 

OF con 
= Dy (9) A) = 0. (60) 
In Example 9, the characteristic projections are given by 


F(z, y; A) 


s+tsin(s) —2 =0, (61) 


where s is constant along each characteristic projection. Their envelope is found by differen- 
tiating with respect to s, 
1+ tcos(s) = 0, (62) 


which occurs when the Jacobian O(x, y)/O(s, 7) is zero. 
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Increasing ¢ 


-1+ 


Figure 6: The solution u(x,t) from Example 9 plotted versus x for t = 0.5, 1.0,1.5,2.0. The 
initial solution u(x,t) = sin(#) is shown as the thicker line. 


3 Weak solutions and shocks for first order PDEs 


3.1 Solutions with discontinuous first derivatives 


In a so-called classical solution, u is smooth so that its first derivatives 0u/Ox and Ou/Oy 
are continuous, and the partial differential equation (5) is satisfied at each point in the (z, y) 
plane where the solution is defined. Now we consider the possibility that the first derivatives 
of u might be discontinuous across some curve C in the (x,y) plane. The idea is to patch 
together classical solutions on either side of C although, on C itself, 0u/Ox and Ou/Oy are 
not well defined. 

Suppose that C' is parameterised by x = 2(&),y = y(&€). We use the superscript * to 
denote the solution on either side of C. By differentiating both solutions along C' we find 


dut Out dx = Out dy du- Ou" dx Ou™_ dy 


dé du d&* Ody dé’ dé Oa d€ Oy dé (63) 


Although the first derivatives of u are discontinuous across C, wu itself is assumed to be 
continuous, so ut = u-. It follows that du*/d€ = du7 /d€é and therefore 


da [u]* dy [Ou] _ 
dé Ean a Fie ef) 


where [f]t = f+ — f~ represents the jump across C. 
Since u~ are both classical solutions of the partial differential equation (5), we have 


Ou* Ou~ 
a~ — +b*¥—=c. 65 
Ox Oy a6) 
Recall that u is continuous across C and, therefore, so are a, b and c: at = a~ =a and so 


forth. By subtracting the equations on either side of c, we thus find 


[4] +0[%] =o (66) 
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In (64) and (66), we have a homogeneous linear system for [Ou/Oz]* and [Ou/Oy]|*, which 
must therefore be zero unless the determinant of the system is zero. In other words, the first 
derivatives must be continuous unless 


b— —a— =0. (67) 


This is identical to the equation for a characteristic projection. Thus, the first derivatives of 
u may only be discontinuous across a characteristic projection. Indeed, this may be used as 
an alternative definition of what a characteristic projection is: it is a curve across which the 
first derivatives of u may be discontinuous. 


Example 10 Consider the partial differential equation 


Ou Ou 
| 68 
cog (68) 
subject to the boundary condition 
0 «<0 
u(x,0) = , 69 
(7,0) ‘° a (69) 
The characteristic equations 
dx = dy = du (70) 
give the general solution 
u=xtf(e—y). (71) 
The boundary condition gives 
—s s<0 
_ y 72 
fis) . a (72) 
and the solution is therefore 
cr L£>y 


Notice that the first derivatives of u are discontinuous across the characteristic y = x that passes 
through the origin, but u itself is continuous. 


3.2 Weak solutions 


In the previous section we showed that classical solutions may be patched together in such a 
way that the first derivatives of u are discontinuous across a characteristic projection. Now 
we attempt to do the same for solutions in which wu itself has a jump across some curve 
in the (x,y) plane. Selecting a unique solution is inherently more problematic when wu is 
discontinuous, as the following example illustrates. 


Example 11 Consider the partial differential equation 


Ou Ou 
a + Dy = 0, (74) 


subject to 


0 y<0, 
won ={' ae (75) 
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Figure 7: Schematics of the possible characteristic projections near a discontinuity in u at 
the point s = sg on T. 


Now we try to find a curve y = f(x) such that the solution is 
(76) 


This clearly satisfies the partial differential equation, everywhere except on y = f(a), and the initial 
condition so long as f(0) =0. Otherwise there does not appear to be any unique way of choosing f. 
Note, though, that u is constant along each of the characteristic projections which, for this linear 
partial differential equation, are given by y= x-+ const. We therefore have u = 0 on the characteristic 
projections leaving x = 0,y < 0, and u = 1 on those that come from x = 0, y > 0. This implies that 
u=Oiny<a2andu=1iny>4a2, +e. that the correct choice is f = x. 


Example 11 illustrates a plausible way of selecting a unique solution for semilinear equa- 
tions, for which the characteristic projections are determined independently of the solution. 
Suppose the initial data have a discontinuity at some point s = sg along I’. On either side 
of s = so, a unique solution is determined on each characteristic projection leaving [. This 
suggests that the discontinuity in u simply propagates along the characteristic projection 
through so (e.g. the line y = x in Example 11). 

Unfortunately this approach does not work for quasilinear equations. The problem is that 
the characteristic projections cannot be found in advance: they depend on the solution. If u 
has a discontinuity at s = 59, then so does the slope of the characteristic projections leaving 
I either side of sg. The two possible outcomes are illustrated schematically in Figure 7. In 
diagram (a), the slopes of the characteristic projections leaving either side of s = sg are such 
that they diverge. There is therefore a region between the limiting characteristic projections 
in which we do not know the solution. In Figure 7(b), the characteristic projections from 
either side of s = sg cross, so there is a region in which they overlap and in which there are 
therefore two possible solutions for uw. 

To resolve these difficulties, we now reformulate the problem in such a way that it makes 
sense even if u is discontinuous. The idea is to turn the partial differential equation (5) 
into an integral equation since, although a discontinuous function does not have well-defined 
derivatives, it may readily be integrated. The first step is to write the partial differential 
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Figure 8: Schematic showing the boundary curve I’, closed by a curve ¥y to enclose a region D. 


equation in so-called conservation form 


dP AQ 
eta (77) 


where P, @ and R are functions of x, y and u. 


Example 12 The semilinear equation 


Ou Ou 
a(z,y) 5 + (2, y)a- = ez, 4,4) (78) 
at y 
may be rewritten as 
O Oa O Ob 
a) usd mies UG =; (79) 


and thus in conservation form (77), with P = au, Q = bu, R=c+uda/Ox + udb/Ody. 
Example 13 The inviscid Burgers equation 

+u— =0 (80) 
may be written in conservation form with P =u, Q = u?/2, R=0. 


Suppose that u is given as usual on some curve [ in the (x,y) plane, and we wish to 
determine the solution for u in some domain D, formed by closing I with an additional curve 7, 
as illustrated in Figure 8. Now, we multiply (77) through by a continuously differentiable test 
function w, assumed to vanish on y, to obtain 


OP OQ. 
On a a = Ry, (81) 
which may be rewritten in the form 
O 0 _ p0v Ow 
ag PY) + (QU) = PZ +5" + Ry. (82) 


Now we integrate both sides of this equation over the region D: 


a) a _ Ov HOW 
II ag PY) + py (Qe) dady = I ees 25, + Ry dady. (83) 
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Figure 9: Schematic showing the shock C’ dividing D into two regions D, and Dy. The 
integration paths on either side of C are denoted Cy and Co. 


We apply Green’s theorem to the left-hand side and use the fact that w~ is assumed to vanish 
on ¥: 


[wv Pdy—Qdz) = II PIE + QE + Reddy. (84) 


A so-called weak solution of the partial differential equation (77) is a function wu that satis- 
fies (84) for any suitably differentiable test function w. If u is continuously differentiable, then 
the steps that led from (77) to (84) may be reversed. Thus, any continuously differentiable 
u satisfying (84) is a classical solution of (77). However, (84) makes sense when u is non- 
differentiable or even discontinuous, while the original partial differential equation (77) does 
not. This is because, by using Green’s theorem, we have removed the need to differentiate u: 
only the test function is differentiated. 


3.3 Shocks 


Now we show how the weak formulation (84) allows us to make sense of solutions in which u 
is discontinuous across some curve C’ in the (x,y) plane. Such a curve is called a shock; this 
name arises from the occurrence of such solutions in gas dynamics. If the shock is initiated 
on I, then it will divide our integration domain D into two regions D, and D2, as shown in 
Figure 9. Thus the area integral on the right-hand side of (84) can be written as 


[P53 Fe $Opy, + Roaedy = ff Po 3, + Rvdedy 
“7 poe Qe ° + Ry dandy, (85) 


Now, inside each of D; and Do, the solution is supposed to be continuously differentiable, 
so we can write 


} a, zt Ru dady 
dP AQ 


= ff serv + FQ) +0 (R- FE - FP) aody, 6) 
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where i = 1 or 2. The term in brackets is identically zero, because of (77), and Green’s 
theorem may be applied to the remainder to give 


IL Poe 5, + Rvdedy =f o(P dy — Qdz). (87) 


As indicated in Figure 9, the integration curves 0D; and 0D2 comprise sections of I and y 
joined to curves C; and C2 adjacent to the shock on either side. When the two integrals are 
summed, the result is 


¢ w(Pdy -Qar) = f w(P dy — Qaa), (38) 
OD\+0D2 T+C,—C2 


since w is zero on y. Notice the difference in sign because C; and C2 are traversed in different 
directions. The right-hand side of (84) may therefore be written in the form 


[f?5 poe +Q5° + Rudedy = f UtPin Oday. (89) 
D [+C1—-C2 
The integral along TI’ cancels with the left-hand side of (84), and we are left with 
p(Pdy—Qdzx)— | p(Pdy- Qdzx) =0, (90) 
Ci C2 
or 
[ viertay - att ae) = 0, (91) 


where [F]* denotes the jump in F across the shock C. 
Since this relation must hold for any (suitably smooth) test function 7, the term in 
brackets must be identically zero and we therefore obtain 
dy — [Q\- 


re Pt (92) 


This so-called Rankine—Hugoniot condition determines the slope of the shock in terms of the 
discontinuities in P and Q across it. 

For semilinear equations, as shown in Example 12, we have P = a(x, y)u, Q = b(a, y)u so 
(92) reduces to 


= = (93) 


This is identical to the slope of characteristic projections so, for semilinear equations, solutions 
may only be discontinuous across characteristic projections. Thus the solution obtained in 
Example 11 is a valid weak solution. For general quasilinear equations, shocks are different 
from characteristic projections. 


Example 14 Find a weak solution of the problem 


Ou Ou 

scfiae mo t 4 

ap Uae 0, > 0, (94) 
1 
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Shock 


Figure 10: The characteristic projections and shock for Example 14. 


We look for a solution in which u is piecewise constant: u = 1 fora < X(t) and u = 0 for 
x > X(t). Our only remaining task is to determine the position X(t) of the shock. In conservation 
form we have 


Ou O 7, a 
Het Bg (2™) = 9 (96) 
so the Rankine-Hugoniot condition (92) gives 
+ 
ax (elt tue o 
d¢° ft #2 2 
Hence X(t) =t/2 and the solution is 
1 a< st 
u= (98) 
0 «> st. 


The characteristic projections have slope given by da/dt = u. They are illustrated in Figure 10; 
note that the characteristic projections travel into the shock from either side. 


The great advantage of allowing weak solutions in which u is discontinuous is that we 
can eliminate multiple-valued solutions like that shown in Figure 6. As soon as the Jacobian 
becomes zero, we insert a shock that prevents the characteristic projections from crossing. In 
principle, the Rankine—Hugoniot condition (92) determines the position of the shock, and we 
can find the solution on either side by the usual characteristic methods. 


3.4 Nonuniqueness of weak solutions 


We have shown how one may construct weak solutions by patching together classical solutions 
and applying the Rankine-Hugoniot condition (92) across any shocks where u is discontinuous. 
This allows us to avoid the unphysical possibility of u becoming a multiple-valued function. 
Unfortunately (92) is in general not enough to determine the solution uniquely. There are 
two further problems to be addressed. 
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Problem 1 


The Rankine-Hugoniot condition (92) depends on the functions P and Q that appear in the 
conservation form (77). However, there may be many different ways of expressing the same 
PDE in conservation form, and a different conservation form leads to a different Rankine— 
Hugoniot condition. 


Example 15 The PDE 


nas avi) 99 
Ox “Oy 29) 
may be written in many different conservation forms, including 
Ou O 2 ) 2 6) 3 
ape or 5g (2) yo (100) 
These lead to the two different Rankine-Hugoniot conditions 
+ salle 
dy _ [du2]* uy tue 7 dy — [zu] 2 (tc. + ULU + u?) (101) 
dx fujJt 2 dx [du2]" 7 3 (uz + u_) , 
respectively. 
Solution 


The key is to make sure that the functions P and Q in the conservation law correspond to 
real physical quantities (e.g. mass, momentum, energy, etc.). Then the Rankine—-Hugoniot 
condition (92) ensures that these are conserved across shocks. 


Example 16 Consider the following model for traffic flow. Let u(x,t) be the density of cars on a 
stretch of road, where x is distance along the road and t is time. Suppose that the speed v of each car 
is related to the local density by v = (1—u). Then the flux of cars is given by uv = u(1— u) and the 
equation representing conservation of cars is 


— + =—(u(1—u)) =0. (102) 


(103) 


fe(@—a-w)] 20 ap 


This equation ensures that the rate at which cars enter a shock from one side equals the rate at which 
they exit the other. 
If (102) is rewritten as 


which may be rearranged to give 


O12 O12 23 

ai (au) + 5, (au — gu") = 0, (105) 
then the corresponding Rankine-Hugoniot condition has no physical interpretation, and the net flux of 
cars would not be preserved across a shock. 
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Problem 2 


The strategy described above selects a particular weak formulation (7.e. a particular choice 
of P and Q) from the many that may be available. The second problem is that this weak 
formulation may still admit many solutions. 


Example 17 The problem 


Ou Ou 0 «£<0 
pd eres = 106 
at “ag em) ‘ r>0 es 
has (at least) two possible weak solutions, namely 
0 z<0 
0 <t/2 
u=fa/t O<aK<t and nat as (107) 
1 x>t 


It is readily verified that 0, «/t and 1 all satisfy the PDE, and the discontinuities in the first derivatives 
of u, occur across characteristic projections, so u, is a valid solution. Similarly, uz clearly satisfies 
the PDE on either side of the shock at x = t/2, and the Rankine-Hugoniot condition da/dt = 5 is 
satisfied, s0 Ug is also a valid weak solution. 


There are several ways around the problem illustrated in Example 17, including the fol- 
lowing. 


1. Entropy 
Here we pose a function of u that must increase across a shock. 


2. Viscosity 
This involves introducing a higher derivative to regularise the equation (i.e. smooth 
out the shock). 


Example 18 As shown in Example 17, the equation 


Ou Ou 
=) + uz =0 (108) 


+u— =e (109) 


has a unique solution, given suitable initial and boundary data, for any positive value of €. So a unique 
weak solution of (108) may be selected by solving (109) and then letting € > 0. 


3. Causality 
This means ensuring that information travels into the shock, not out of it. 


It may be shown that each of these approaches results in the same unique weak solution 
being selected. We only use the third option, which does not require us to bring any more 
physics into the problem, apart from a recognition that one variable (t) represents time. 
Information travels along characteristic projections, starting from the initial data, in the 
direction of increasing t. A shock solution is causal only if this information travels into the 
shock from either side (as in Figures 10 and ??). By disallowing shocks that do not satisfy 
this condition, we narrow down the possible weak solutions to just one. 
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(a) (b) 


AAAAAAAKAKA 


Figure 11: Characteristic projections for the two possible solutions, (a) u; and (b) ug, to 
Example 17. 


Example 19 The characteristic projections for the two possible solution u, and ug to Example 17 
are shown in Figure 11. Notice that the shock solution in Figure 11(b) has characteristic projections 
travelling out of the shock. Thus this solution is not causal and must be discarded. The solution uz 
displayed in Figure 11(a) is the correct weak solution. 


In Example 19, the characteristic projections have slope dr/dt = u, while the slope 
of the shock is given by the Rankine-Hugoniot condition da/dt = 5 (U4. +u_). Thus the 
characteristic projections travel into the shock from either side so long as 


u— > 5(u4 +U_) < us, (110) 


which may be rearranged to give 
. > UL. (111) 


This is the condition for a shock to be causal: u must be greater behind the shock than it is 
ahead. This is why the shock solution in Example 14 is acceptable while that in Example 19 
is not. 
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4 First order nonlinear equations 


In this chapter we introduce a method that will enable us to solve fully nonlinear first order 
PDEs by solving an associated system of ODEs (Charpit’s equations). 

We consider general first-order nonlinear scalar PDEs, i.e. not (necessarily) quasi-linear. 
The general form of such an equation is 


F(p,q,u,2,y) = 0, (112) 
where we use a e. 
dn op q (113) 
as shorthand, so that 
a = = (114) 


The case of quasilinear equations corresponds to F’ being a linear function of p and q, i.e. 
F(p, g,U,2,Y, ) = a(x, Y; u)p + b(a, Y; u)q = ela, Y; u). (115) 
4.1 Charpit’s equations 


If we differentiate (112) with respect to x and y, we obtain 
OF Op OF 0q _ OF OF 


Op Ox =Oqg Ou Ox Pau’ lisa) 
OF Op OF 0q _ OF OF 
Op Oy = Oq Oy —— Oy T Bu? tah) 
or, using (114), 
OF Op OF Op. OF OF 
Op Ox =Oq Oy Ox Pou’ ane 
OF 0q OF 0q OF OF 
Op Ox = Oq Oy — Oy T Oy oe) 
So, if we define characteristics or rays as curves (alr); y(T)) satisfying 
dx OF dy OF 
dr Op’ dr 0q oe 
then, along these curves, 
dp OF OF dq OF aF 
dr On? Ou’ dr Oy 1 Ou uy) 


We therefore have a system of four ODEs for z, y, p and q satisfied along the rays. Recall, 
though, that in general F’ depends on wu also, so to close the system we also need an ODE for 
u along the rays, namely 


du Oudx Oudy OF OF 
= = : 12 
dr Oxdr Oydr - Op +4 Og aa 
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In summary, we have the following system of ODEs for x, y, p, q and u, known as Charpit’s 
equations: 


“ = > (121a) 

wu = 7 (121b) 

be ae ome 

— ~ ae (1214) 

“ =P 5, aay (121e) 
It is easily verified that these reduce to the usual characteristic equations 


for quasi-linear equations where F’ takes the form (115). 


4.2 Boundary data 


As for quasilinear scalar equations, Cauchy data specifies u along some curve [I in the 
(x, y)-plane: 
t= r9(s), y= yo(s), U= uo(s), (123) 


for s in some (possibly infinite) interval. We also require initial conditions for p and q, say 
p = po(s), ¢ = qo(s), which are obtained by differentiating ug with respect to s and using the 
PDE (112): 


= | .  F(po,40,%t0,.£0, Yo) = 0. 124 
qs 7 Pog, + 0G, (Po; Go, Uo, Lo, Yo) (124) 


By the implicit function theorem, the two equations (124) may be solved uniquely (in principle, 
if not explicitly) for pp and go provided the condition 


dxo OF dyo OF 
ds Oqg ds Opo 


is satisfied. This is the same as insisting that I not be parallel to a ray. 

Charpit’s method consists of solving the ODEs (121) for (p,q,u,x,y), with (123) and 
(124) as initial data at 7 = 0. This gives (p,q,u,2,y) all as functions of s and 7 and, in 
principle, allows us to reconstruct the solution surface u = u(z, y). 


(125) 


Example 20 Sugar on a spoon 

Consider sugar piled up on a spoon such that its height is given by u(x,y). At criticality, just before 
the sugar would start to slide off the spoon, the angle between the normal to the surface and the vertical 
(0,0, 1) ts a prescribed constant, y, the angle of friction. That is, 


2 2 
0 
(=) + (55) +1=sec? (126) 
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After rearranging and normalising, this can be written as the Eikonal equation 


du\?  (du\* 
— —]}] =1 127 
(a) + (a) =» an 
which is of the form (112) with 
F(p,q) = 3 (e’ +¢° - 1). (128) 
Charpit’s equations for this particular F' are 
dx dy dp dq du 2, 2 
dr dr” dr = dr ” dr” +4 2) 


Notice that p and q are constant along rays and, hence, given by their boundary values: 
P = po(s), q = 40(8). (130) 
The remaining ODEs are then readily integrated to give 
« = x0(s) + po(s)r, y = yo(s) + go(s)r, u = uo(s) +7. (131) 


Notice that the slope of a ray is given by qo(s)/po(s) which is constant along each ray. Thus the rays 
are straight lines, along which u increases linearly with T. 
At the edge of the spoon, the height is zero, so up(s) =0. Then we have the system 


dzo dyo D) 2 
qs 20+ Gy 0 =O Po + % (132) 
for po and qo, whose solution is 
ty! xt 
Po = at go = 0 (133) 


(26)? + (46)? (x)? + (yo)? 


where’ is used as shorthand for d/ds. The vector (po, qo) is the unit normal to the boundary T. Hence 
the rays are straight lines perpendicular to T and u(x,y) is simply the distance of the point (x,y) 
from T. 

Notice that there are two possible solutions corresponding to the + in (133). The correct solution 
is chosen by ensuring that the rays propagate into the region of interest, not out of it. So, in (133), 
we have to choose (po, qo) to be the inward pointing normal. Otherwise the solution corresponds to the 
sandpile outside a spoon-shaped hole in a table. 


4.3. Proof that Charpit’s method works 
If we differentiate F along a ray, we find that 


dF OF dp OF dq | OF du OF dx , OF dy _ 
dr Opdr  Oqdr Oudr = Oxdr Oydr — 


0. (134) 


Since the boundary condition (124) sets F to zero on the initial curve T’, it must therefore 
be zero everywhere along the rays passing through I’. Hence p, q and u satisfy the equation 
F(p,q,u,£,y) = 0 everywhere in the domain of definition where there are rays emanating 
from [. 

This is not quite sufficient to prove that u is a solution of the original nonlinear PDE. We 
still have to show that the functions p and q that result from solving Charpit’s equations are 
equal to 0u/Ox and Ou/Oy respectively. To do this, we first prove that ¢ = 0, where 


Ou Ox Oy 
ce: 


(135) 
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(b) 


Figure 12: (a) Rays for a sugar heap on an elliptical spoon with a = 2 and b = 1; the bold 
line marks the ridge. (b) The corresponding pile height u(z, y). 


By differentiating ¢ with respect to 7 and rearranging, we obtain 


Ob OF OF oO (> OF OF ) 


a pe ae Op ae (136) 


du Os 


The final term is identically zero by (121), and we have already shown that F = 0 in the 
domain of definition, which implies that OF'/0s = 0. Hence ¢ satisfies 


=o (137) 


with, by (124), 6 =0 at r = 0. Provided OF'/Ou is bounded, it follows from Picard’s theorem 
that ¢ = 0 in the domain of definition. 
From this fact and Charpit’s equation for u, we obtain two simultaneous equations for 


Ou/Ox and Ou/Oy: 
Ox Oy Ou  OxOu  Oydu 


ce a Or! Or Orda OF Oy’ vey 
x Oy Ou OxOu_ Oy Ou 
as? + as’ Os  0sdxn Os Oy” tsa) 
This system has a unique solution, namely 
Ou Ou 
———s = —_ 1 


provided the determinant of the right-hand side is nonzero, and this determinant is the Ja- 
cobian of the transformation from (s,7) to (x,y), that is 
_ Oyox OxOdy Ox0F OyOF 
070s Or0s Os 0q Os Op’ 
Hence we have shown that u(x, y) satisfies the nonlinear PDE (112) as long as J is nonzero. 
This proof applies in particular to the special case where F is given by (115), and hence 
establishes that the method of characteristics works for scalar quasi-linear PDEs. 


(140) 
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4.4 Discontinuities 


Example 20 is an example in which the rays intersect. This happens where the Jacobian J 
defined by (140) first becomes zero. Note that this reproduces the criterion (125) for Cauchy 
data not to determine a unique solution on I’. If rays are allowed to cross, then the solution 
becomes multi-valued, which is clearly unphysical for a pile of sugar. Instead, we must allow 
shocks to form across which the solution is discontinuous. Recall that, for nonlinear PDEs, 
shocks are different from characteristics. The conditions that must be applied across a shock 
depend on the physical situation being modelled. For the sugar heap problem, it is clear 
that u must be continuous everywhere, since a discontinuity in u (corresponding to a vertical 
“cliff” ) cannot be sustained. This forces the shock, i.e. the ridge line, to be along the x-axis 
as shown in Figure 12. In general, the region of validity of the solution obtained by Charpit’s 
method is bounded by curves on which J = 0. 


4.5 Geometrical optics 


The propagation of sound or light waves in two spatial dimensions is governed by the wave 
equation 
Or Or aa om) 
Ox2 © Oy? c2 At?’ 
where w is some state variable such as pressure or electric field and c is the wave speed. We 
look for time-periodic (or “monochromatic” ) solutions with constant frequency w by setting 


W(x, y,t) = o(x, ye". (142) 


Then ¢ satisfies the Helmholtz equation 


(141) 


Pb ae 
V7b+ kh? = 5 +.5+k =), 143 
where k = w/c is the wavenumber (i.e. 27 divided by the wavelength). 
The theory of geometrical optics arises from the limit k — oo, which is valid over length- 
scales much longer than the wavelength. To analyse the behaviour of (143) in this limit, we 
use the so-called WKBJ method, which involves writing ¢ in the form 


$(z,y) = A(x, ye), (144) 


where A and u represent the amplitude and phase respectively of the solution. Then (143) 
becomes 


V?A + ik (AV2u + 2VA: Vu) + kA (1 — |Vul?) =0. (145) 
We seek solutions in which A is an asymptotic expansion of the form 
A A 
AAG ss (146) 


At leading order, (145) implies that wu satisfies the Eikonal equation 


|Vul? = 1. (147) 


B5.2 Applied Partial Differential Equations 34 


reflected 


incident 


Figure 13: Illustration of Snell’s law. 


Then the successive terms in the amplitude expansion satisfy the transport equations 
2Vu-VAo + AoV7u = 0, (148a) 
2Vu-VAn + AnV2u = iV7An—1, n>1. (148b) 


The Eikonal equation (147) may be solved exactly as in Example 20. The rays correspond 
to light rays and all the familiar properties of geometrical optics, for example that light travels 
in straight lines, follow from the solution of Charpit’s equations. 


Example 21 Reflecting plane waves 

One obvious solution of (147) is u= 2x, which corresponds to a plane wave moving in the x-direction. 
Now we examine what happens if such a wave impinges on a reflecting wall given by a curve T in the 
(x, y)-plane. We decompose the state variable } into an incident wave ¢;, namely a plane wave with 
constant amplitude a, and a reflected wave gr: 


b= ¢1+ or, orp=te™, (149) 
Now we apply the WKBJ ansatz to or: 
br= Acikuley). (150) 


The boundary conditions depend on the physical situation being modelled and exactly what the state 
variable ¢ represents. The simplest case is to impose the Dirichlet condition 6 = 0 onT, which leads 
to 


Uu=f, Ao = —a onl. (151) 


Other possibilities are 06/On = 0 or the “Robin condition” 0¢/On + A> = 0, but it ts readily verified 
that the leading-order boundary conditions (151) are unchanged in either case. 
The solution of Charpit’s equations for the Eikonal equation was already obtained in Example 20: 


p=pol(s), G=4q0(s), ®t=20(s)+po(s)T, Yy=yo(s)+q0o(s)T, w=uo(s)t+T. (152) 


For simplicity we suppose that s parametrises arc-length along T so we can write x = cos, yj = sind 
where @ is the angle between T and the x-axis. The boundary condition (151) implies that uo(s) = xo(s), 
and then po and qo are obtained from 


LyPo + Yodo = Xo, po +45 = 1. (153) 


This system has two solutions, one of which is po = 1, qo = 0, corresponding to the incident wave. 
The reflected wave is given by the other solution 


po = 1 — 2(y))? = cos(26), Go = 2xpy = sin(26). (154) 


Hence the reflected ray makes an angle of 20 with the x-axis. This is Snell’s law: as illustrated in 
Figure 13, it implies that the angle of incidence to the wall equals the angle of reflection from it. 
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caustic a 


Figure 14: (a) Reflected rays for Example 22. (b) Rays truncated by a caustic on which 
J =0. 


Example 22 The caustic in a teacup 
As a special case of Example 21, we now consider the case where T is the unit circle, parametrised 
(say) by xo(s) = cos(s), yo(s) = sin(s). 

As shown in Figure 14(a), the reflected waves start to cross a finite distance from the circle. The 
envelope of the rays, determined either via the Jacobian condition 


Ox OF OyOF 
= -  — 1 
Os Oq = Os Op u a) 
or by solving the envelope equations 
OF 
F :s) = = 1 
(t,y38) = ay =0, (156) 


defines a curve of concentrated light, as can be observed by shining a light in a teacup. 


A single-valued ray solution may be obtained by truncating rays at any caustic where 
the Jacobian is zero. It may be shown that the asymptotic ansatz (144) breaks down, with 
A — oo as the caustic is approached. The method of matched asymptotic expansions yields 
the appropriate correction in the neighbourhood of a caustic and allows the behaviour in the 
dark zone beyond the caustic (corresponding to complex rays) to be found. 
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5 First order systems 


5.1 Introduction 
In this chapter we consider first-order systems of PDEs of the form 


Ou 


A(x, y,u) Aa 


+ Bley.u) 5 = e(e.y.u), (157) 
where now the unknown wu is an n-dimensional vector function of x and y, c is also an 
n-dimensional vector, A and B are n X n matrices. Again, we only consider quasilinear 
equations, so that (157) is linear in the first derivatives of u. Equations of the form (157) 
arise frequently in physical problems (e.g. gas dynamics), in which it is often the case that 
one of the independent variables represents time, so we will also write (157) in the form 


a) a) 
A(c,t, uw) 4 Be, t,u) 5 = (x,t, u). (158) 
The other motivation for studying (157) is that many higher-order scalar equations may be 
written in this form. 


Example 23 The shallow-water equations describe the flow of a thin layer of liquid over a flat 
surface under gravity. If x measures horizontal distance and t is time, then the dimensionless height 
h(a,t) and velocity u(x,t) satisfy 

Oh Oh Ou Ou f wot x Oh _ 
Ot Ox Ox 


0, (159) 


which may be written in the form (158), with 


ee 


Example 24 Laplace’s equation 


asta =0 (161) 


may be written as a first-order system by setting u = 0¢/0x, v = —O6/Oy. Then u and v satisfy the 
Cauchy—Riemann equations 


Ou Ov Ou Ov 
oie — 4+ = 162 
Ox Oy Ze Oy 2 Ox 0 an02) 


which are in the form of (157), with 


(J. G9 GB Gm 


(Once u and v are known, ¢ is found by solving either 0¢/0x = u or 0¢/Oy = —v.) 


5.2 Cauchy data 
Definition 


For the PDE system (157), Cauchy data is to specify u on a curve [ in the (a, y) plane, i.e. 


xL = x0(s), y = yo(s), u = Ud(s), 81 <8 < 89, (164) 
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and the PDE (157) together with the data (164) is known as the Cauchy problem. As for the 
scalar case, we now ask whether the Cauchy problem determines the first derivatives of wu. 
Differentiating (164) with respect to s, we find 


duo _ dxo Ou 5 dyo Ou 
ds ds Ox ds Oy’ 


Now we consider (157) and (165) as a (2n) x (2n) matrix system for the (2n)-dimensional 


vector (Ou/dx,0u/dy)?: 
A B Ou/dOx\ [ec 
(chr) (Buran) = (2a) = 


where ’ is shorthand for d/ds and I is the n x n identity matrix. The system (166) may be 
solved uniquely for the first derivatives of w provided the determinant of the matrix on the 
left-hand side is nonzero, and this determinant may be rearranged to give 


det (xp B — yy A) #0. (167) 


(165) 


This is the condition on the initial data for the the first derivatives of wu to be locally deter- 
mined. It clearly reduces to the condition found for scalar equations when n = 1. 


Cauchy—Kowalevski theorem 


Now we state a generalisation of the theorem previously introduced for scalar PDEs. For 
simplicity we suppose that a coordinate transformation is used to shift the boundary [I onto 
the y-axis, where we specify u: 


u = u(y) on r=0, ySy<ye. (168) 


Clearly, so long as uo is differentiable, we can calculate Ou/Oy directly: 
Ou duo 


= =0 <y < yo. 169 
Oy dy ee 7 WAYS Ye (169) 
We can then use the PDE (157) to solve for 0u/Oz, 
Ou Ou Ou 
—=A'c-~A'B— = — 170 
Or c Oy f («. Y; U, Oy ) rh Say, ( ) 


so long as A is invertible. Thus the PDE may be written in the form (170) provided |A| 4 0, 
which is the same as condition (167), with ro = 0, yo = s. 

Now suppose that uo(y) is analytic at a point yo € (yi, y2) and that f is analytic in all 
its arguments at the point 


(0.10, uot), on) 


Then the Cauchy-Kowalevski theorem says that the Cauchy problem (168, 170) has a unique 
analytic solution u(x, y) in a neighbourhood of (0, yo). 

The proof of this theorem works by constructing a Taylor expansion for wu and showing 
it has a finite radius of convergence. It may readily be extended to the case of a general 
analytic initial curve [ if the condition (167) is satisfied. However, the result is entirely local: 
the radius of convergence may be extremely small, meaning that wu may become singular or 
nonunique a short distance from the initial data. Moreover, the hypotheses of the theorem 
are rather restrictive: it says nothing about existence or uniqueness of solutions when the 
initial data (say) are nonanalytic. 
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5.3. Characteristics 
Definition 


There is not a straightforward geometric interpretation of the PDE system (157) as there is 
in the scalar case. At the very least we have to imagine n solution surfaces — one for each 
component of wu. There is also no obvious generalisation of the characteristic ODEs. Recall, 
though, that for scalar PDEs, characteristics are curves on which Cauchy data does not allow 
the the first derivatives of u to be determined uniquely, and this definition does generalise to 
n dimensions. 

Consider a curve in the (x,y) plane, parameterised by 7, and suppose that w is known on 
this curve. Then, according to (167), the first derivatives of w may be found locally unless 


daz dy _ 


A curve (x(7T), y(T)) satisfying (171) is called a characteristic projection. Since there is no 
generalisation of the concept of characteristic introduced previously for scalar PDEs, it is 
usual to drop the word “projection” and refer to these plane curves simply as characteristics. 

As for scalar PDEs, an alternative definition is that characteristics are curves across which 
the first derivatives of w may be discontinuous. Suppose wu is continuous across a curve C’ in 
the (x,y) plane, but its first derivatives take different limiting values (denoted by * and ~) 
on either side of C. If C is parametrised by x = x(€), y = y(&). Then, by differentiating wu 
along either side of C, we obtain 


du Out dx Out dy Ou dx | du" dy 
dé Ox dé Oy dé Ox dE dy dé 
and by subtracting these two equations we find 
dz [Ou]* dy [du]* 
= 0. 173 
dé ae |_ +a ae. 4 
Similarly, because u satisfies the PDE (157) on either side of C’, we have 


a [4] +3 (24]" =o an 


(172) 


The homogeneous system (173, 174) implies that the jumps in 0u/Ox and Ou/Oy must be 
zero unless the determinant of the system is zero, which implies that 


d d 
det (Gs = sta) =0, (175) 
i.e. that C is a characteristic. 


Classification 
The slopes of the characteristics satisfy the eigenvalue problem 
= where det (B— A) =0. (176) 


Thus X satisfies an nth-order polynomial equation, whose roots may be complex in general. 
A 2 x 2 system may be classified as follows, depending on the eigenvalues A. 
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e If there are two distinct real eigenvalues, then the system is said to be hyperbolic. 
e If there is one repeated real eigenvalue, then the system is parabolic. 
e If the eigenvalues are complex, then the system is called elliptic. 


Example 25 Consider the quasilinear second-order PDE 


ae, 42 Po | Fb _ 


iL 
“Ax? Oxdy ses Oy? f ee) 


where a, b and f are in general functions of x, y, ¢, 0¢6/0x and 0¢/Oy. We can write (177) as the 
first-order system 


a b\0Ou b c\ Ou (f _ (06/0x 
( On + 6 4 ay (i) ; where w= Gs ; (178) 
Ifa, b, c and f are independent of d, then we can ignore the uncoupled equation 0¢/0x = u, and the 
characteristic slopes satisfy 


eas c—Ab 


1 d [=o 


=> ad—2br’\+c=0. (179) 


The system is thus hyperbolic if b? > ac, parabolic if b? = ac, elliptic if b? < ac. 


In dimensions higher than two, there are clearly many possible combinations of real, 
complex, distinct and repeated roots of the polynomial equation (176), and there is no such 
simple classification. However, we still define an equation as hyperbolic if (176) has n 
distinct real roots A. Since the matrices A and B depend on z, y and, in general, also on the 
solution u, the type of the equation (7.e., hyperbolic, elliptic, parabolic or some hybrid) may 
also vary with position. 

Now, according to the Cauchy—Kowalevski theorem, provided all our coefficients and initial 
data are analytic and the condition (167) is satisfied, there is a unique solution for wu in a 
neighbourhood of [. Nevertheless, unless the PDE is hyperbolic, the Cauchy problem is in 
general ill posed. This may manifest itself in several ways. For example, the unique local 
solution may blow up arbitrarily close to or may be pathalogically sensitive to the initial 
data. 


Example 26 For the Cauchy—Riemann equations, the characteristic slopes satisfy 


-r’ -1 F 
| 1 0 => A=, (180) 
so the system is elliptic. 

Suppose we try to solve (162) in y > 0 subject to the Cauchy data 


6? 


u = uo(x) = 0, v= vo(t) = a 


on y=0, (181) 
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where € and 6 are positive constants. Notice that the initial data are analytic and bounded: vo(ax) is in 
the range (0,€] for all x. It may be shown! that the solution of this problem is 

_ 257. ry = <6? [a? — y? + 6°] 

fa? + (y — 6)7] fa? + y+ 6)9)’ "te + (y — O22 + (YF 5] 


Thus, however small € is, both u and v blow up at the point (0,6): by choosing 6 small, we can make 
the solution break down arbitrarily close to the initial curve y = 0. 


(182) 


For the remainder of this chapter we restrict our attention to hyperbolic systems, for which 
the Cauchy problem is generally well posed, and for which characteristic methods analogous 
to those used for scalar equations can be applied. So, at each point in the (z,y) plane, we 
assume that (171) defines n distinct real eigenvalues A. Thus, by solving dy/da = X for each 
of these n characteristic slopes, we can obtain in principle n families of characteristics for an 
n-dimensional hyperbolic system. 


5.4 Integration along characteristics 


Reduction to an ODE 


Suppose A is a real eigenvalue of (176); recall that is in general a function of x, y and u, 
since A and B are. Now the matrix (B— 4) is singular, so there exists a left eigenvector I", 
such that 


I?(B— XA) =07, that is 7B =X" A. (183) 
Multiplying the PDE (157) on the left by I”, we obtain 
Ou Ou 
LA 2 B= 
Ox 7 Oy . 
Ou Ou 
Al =) 1c 184 
as ( Ox = Oy ) 7 ae) 
Along characteristics, whose slope is dy/dx = A, we have 
du 
PA =e 185 
In c (185) 


This is the equivalent of the ODE satisfied by u along characteristics in the scalar case. There 
is one ODE of the form (185) satisfied along each of the n families of characteristics. 


Riemann invariants 


Since (185) is a single differential equation for the n components of wu, it is usually not 
integrable. However, there are special cases in which the ODE (185) may be rearranged to 
take the form 


< [R(x,y,u)] = 0. (186) 


‘the Cauchy—Riemann equations imply that u + iv is a function of z = «+ iy; here, wu and v are the real 
and imaginary parts of the complex function 


667i 


22+ 62° 


fa= 
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If so, the function R is conserved along the characteristics satisfying dy/da = 4, and is called 
a Riemann invariant. 


Example 27 Consider the system 
which may be written in the form (157) with 


OC 


The characteristic slopes satisfy 


fe a =0 => A=41, (189) 


so the system is hyperbolic. Now we have to solve 


I7(B—\A)=17 G |) = 07, (190) 


and a suitable left eigenvector is I" = (1, +1). When we multiply the system 
ie 


GCC em 


on the left by I”, we obtain 


0 fu O (u 
1 +1) — +1 -1) — = 0. 192 
0 FD ae (s) +e 5, (%) ie 
Using the fact that dy/daz = \ = +1, this may be rearranged to give 
d 
—(uFv) =0, (193) 
dx 
so the Riemann invariants u + v are conserved along the characteristics dy/da = +1. 


In this simple case, the Riemann invariants may be used to write down the general solution. Since 
utv is constant when y+ is constant, we must have u+v = f(y+a2) for some function f. Similarly, 
u—v must take the form g(y — x), so the general solution is 


u=Af(yt+a2)+4g(y—2), v=5f(yt+2x)—Z9(y-2), (194) 


where f and g are two arbitrary functions. Note that, if we cross-differentiate to eliminate v, we find 
that u satisfies the wave equation 


—_=0. (195) 


For linear hyperbolic PDEs with c = 0, as in Example 27, we can always find a com- 
plete set of mn Riemann invariants. Furthermore, for linear PDEs, the characteristics may 
be found independently of the solution. We thus obtain a system of n algebraic equations 
for the components of wu in terms of arbitrary functions that are constant along each family 
of characteristics. This suggests a plausible method for solving hyperbolic systems numeri- 
cally. If A, B and ec are approximated as being locally constant near [’, then the resulting 
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autonomous linear system has a complete set of Riemann invariants. Thus the solution u a 
small distance from I may be found by solving the resulting system of algebraic equations. 
By repeating this process, the solution may be continued further still from the initial data. 
This proposed procedure suggests that the Cauchy problem should usually be well posed for 
hyperbolic systems. 

The following two examples show that, when c is nonzero, even linear PDEs have no 
Riemann invariants in general. 


Example 28 The system 


Ou Ov Ou Ov 
—-> = —- ;=0 196 
Ox Oy ee Oy Ox me 
may be written as 
1 0 )\ du 4 0 -1\0u_ (uty 
0 -1/) Ox 1 0 Oy ~ 0 ; 
where u = (u,v)?. As in Example 27, the characteristic directions are dy/dx = \ = +1 and the 
corresponding left eigenvectors are I? = (1,41). Now, the ODEs satisfied along the characteristics are 
d d 
re (u-v)=utv on iL = El; (197) 
only one of which is integrable: 
d d 
—(utv)=ut+v => e *(u+v) = const Ghee 22 (198) 
dx dx 


Nevertheless, we can use this single Riemann invariant to find the general solution in this case. 
Since e~*(u+v) is constant when (y+) is constant, we may write 


utv=e"f(y+2), (199) 


where f is an arbitrary function. Along the other family of characteristics we have (y — x) = const = k, 
Say, SO 


d 
qg eT US utvaeP fy ta) =e flk + 2x) on y-a=k. (200) 
x 
This may be integrated with respect to x to give 
1 af 2s 
U-v= Te(e—y)/ | e°/? f(s) ds + g(y — 2), (201) 
0 
where g is a second arbitrary function. The general solution is, therefore, 
Tg 1 Loewe f° os 
u=-e"f(yt+2)+=g9(y-—x2)+—e 4 e*/* f(s) ds, (202a) 
2 2 4 
1 1 L ane | a6 
v= xe" f(y+a2)— =g(y—2)—-e” 4 e*/* f(s) ds. (202b) 
2 2 4 , 
Example 29 The system 
Ou Ov Ou Ov 
—- = —- —=0 203 
Ox Oy = Oy Ox ~~) 
has characteristic equations 
d d 
re (u-v) =u on 7s = El, (204) 


which cannot be integrated. There are no Riemann invariants and no way to find the general solution 
explicitly. 


B5.2 Applied Partial Differential Equations 43 


The situation is even worse for fully nonlinear systems, where the characteristics depend 
on the solution wu. Even when such systems have a complete set of n Riemann invariants, 
since we do not know in advance the curves along which each is conserved, we cannot in 
general find an explicit solution. 


Example 30 The shallow-water equations (159) have characteristic slopes given by 


det (B- A) =|"7* yi /atunaP nao 
da 
—==utvh 2 
qT A=t Vh, (205) 


and the corresponding left eigenvectors are (1, +WVh). The characteristic ODEs are 


dh du dz 
vig = on an utvh 
=> £ (u£2Vh) =0 in ale (206) 
dt dé 


so the Riemann invariants ut+2WVh are preserved along the characteristics dx/dt = u+WVh. Although 
the system has two Riemann invariants, we cannot infer the general solution. 


Regions of influence 


Recall that, for scalar PDEs, where Cauchy data are only given on a finite initial curve I’, the 
solution is only determined in the so-called domain of definition, penetrated by characteristics 
emanating from I’. In n dimensions, the domain of definition is the region penetrated by all 
n families of characteristics originating at [. Where there are at least one but fewer than n 
families of characteristics, the solution is influenced but not fully determined by the Cauchy 
data on I. The region swept out by all the characteristics intersecting [ is therefore called 
the region of influence. 


Example 31 We return to the system considered in Example 27, which has the general solution 


u=4f(yta)t+4g(y—-2), v= 5f(yt+2)— Z9(y—2). (207) 


If we are given the Cauchy data u = uo(x), v = vo(x) ony =0,0< a <1, then we can determine the 
two arbitrary functions f and g: 


f(s) = uo(s) + u0(s), O<s<1, (208a) 
g(s) = uo(—s) — vo(—8), l<s<0. (208b) 

The solution for u is, therefore, 
u= z{uoly +x) + voy + x) + uo(% — y) — vo(a—y)f, (209) 


and a similar expression for v may readily be found. 

In (209), the first two terms correspond to the function f and are defined in0 < y+a <1, 
while the final two terms, corresponding to g, are defined inO0<a—y< 1. The region of influence 
is the union of these two domains, where either f or g is defined. The domain of definition is the 
intersection, where both f and g are defined, that is, 0 < y < max(x,1— 2x). They are both illustrated 
in Figure 15. 
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Region of 
influence 


Domain of definition 


Figure 15: Region of influence and domain of definition for Example 31. 


Simple waves 


For nonlinear systems, even if a complete set of Riemann invariants exists, it is usually not 
possible to construct an explicit solution. However, some progress can be made if the boundary 
conditions are such that one Riemann invariant is constant everywhere (7.e. not just along 
characteristics). If this occurs, then we have a functional relation R(x, y,u) = const, which 
allows us to eliminate one of the components of wu and thus reduce the dimension of the 
system by one. For example, two-dimensional systems become scalar PDEs, which may then 
be solved using the methods from chapter 1. These special solutions are known as simple 
waves. 


Example 32 We return to the shallow-water equations (159). As shown in Example 30, these have 
Riemann invariants u42Vh. Suppose we are given u = uo(a), h = uo(x)?/4 ont = 0. Then u—2Vh 
is zero ont =0 and, since it is preserved along the characteristics da/dt = u — Vh, must therefore be 
zero everywhere (in the domain of definition). So we can substitute h = u?/4 into (159) to obtain 


Ou 3udu 


ee 210 
dt’ 2 On . iy) 

which is readily solved to give the implicit solution 
U = Uo(x — 3ut/2). (211) 


5.5 Weak solutions 
Formulation 


As for scalar PDEs, nonlinear hyperbolic systems may have solutions that lose uniqueness 
a finite distance from the initial data. This is illustrated by the solution of Example 32, 
for which 0u/Ox becomes unbounded in finite time if up is ever negative. To continue such 
solutions, it is necessary to allow u to be discontinuous across curves in the (x, y) plane, again 
referred to as shocks. Since the PDE (157) does not make sense on such a curve, we have to 
use a weak formulation of the problem. The theory is very similar to the scalar case, so we 
omit most of the details. 
The first step is to write the system in conservation form 


a eR, (212) 
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Figure 16: Schematic showing the boundary curve I’, closed by a curve ¥ to enclose a region 
D. 


Figure 17: Schematic showing the shock C' dividing D into two regions D; and Dj. The 
integration paths on either side of C are denoted C; and Co. 


where P, Q and R are vector-valued functions of x, y and u.? Now, as illustrated in Figure 16, 
we form a closed region D by closing I’ with a second curve 7, then multiply (212) through 
by a test function ~, assumed to be suitably differentiable and to vanish on y. Then we 
integrate over D and, just as for the scalar case, Green’s theorem leads to the following weak 
formulation of (157): 


_ Ow Ow 
[vay az) = ff PSE + OF + Ro dady. (213) 


A function u(z,y) that satisfies (213) for all suitable test functions w is called a weak 
solution of (157). If w is continuously differentiable and satisfies (213), then it is also a 
classical solution of (157). However, (213) also makes sense if u is discontinuous. 


Shocks 


Now we look for a weak solution in which u is smooth everywhere except a curve C’, across 
which it is discontinuous. As shown in Figure 17, C' divides the region D into two sub-regions 
D, and D2. The integral on the right-hand side of (213) may be split up into two integrals 
over D; and D2 respectively. Since u is smooth within D,; and within Dz, Green’s theorem 


?In fact, an arbitrary PDE system cannot always be written in this form, but it is usually possible for 
physically-motivated problems that are based on conservation laws. 
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may then be used, along with the fact that wu satisfies the PDE (157), giving 


i PS QS + Rvdrdy = us PS QT + Rvardy 


+ ff P Pp +a « mya 


= wv (P dy — Qdz)+ vw (P dy —Qdz). (214) 
OD, OD2 


Then, since w is assumed to be zero on y, (213) reduces to 
[oP ay — Qi? ax) =o. (215) 


where [ |* denotes the jump across the shock. This holds for all test functions 7, and the 
slope of the shock must therefore satisfy the Rankine—Hugoniot condition 


St _ git, (216) 


The scalar Rankine-Hugoniot condition is clearly reproduced if n = 1 but, in higher 
dimensions, (216) gives us n relations between dy/dz and the jumps in the n components of 
u. For semilinear equations, we have 


P= Au, Q = Bu, R=c+ u- U, (217) 
Ox Oy 
so the Rankine—Hugoniot condition is 
dy dy 
+09 _ + _ oy iss 21 
[Au]" re [Bu]" => (2 in 4) [u]" =O (218) 


Thus wu can only be discontinuous if the determinant of the matrix on the left-hand side is 
zero, which implies that dy/dz is equal to a characteristic slope. In other words, shocks occur 
on characteristics for semilinear equations. This is not true for general quasilinear systems, 
though. 


Example 33 The shallow-water equations (159) may be written in the conservation form 


Oh O 7) ) 
ar tan (hu) = = (hu) + = (hu? + $h?) = 21 
at og =e og ae ee 
: OP OQ h hu 
“art Oe =" where P= (1), @= (etait): ee 
so the Rankine-Hugoniot condition is 
+ + 
" ie = 4 2 : (220) 
hu = dt hu + xh = 


If the speed of the shock is denoted by v = da/dt, then this may be rearranged to give 


[(u— vA] * =0, [h(u— v)? + 3h?) =0, (221) 
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Figure 18: Schematic of a bore (with u4 = 0). 


which represent conservation of mass and momentum respectively across the shock. 

A shock solution of the shallow-water equations may be used to model a bore in a river, as depicted 
in Figure 18. Suppose the fluz q coming down the bore and the water height h4 ahead of the bore are 
given, and that the water is stationary ahead of the bore so that us =0. Then we have 


u_-h_ =4, (u_ —v)h_ =—vh,, h_(u-—v)? + hh? =hyv? + hi, (222) 


which constitute three equations to determine h_, u_ (the height and velocity behind the bore) and v 
(the speed at which the bore propagates). 


As for scalar PDEs, there may be many possible ways of writing a quasilinear system in 
conservation form. To obtain a sensible weak formulation of the problem, one should choose 
a form such that P and Q represent physical quantities that one expects to be conserved 
across the shock. Even following this principle, there may be several alternative formulations 
corresponding to different physical scenarios. 


Example 34 The shallow-water equations may also be written in the alternative conservation form 


ae o (dhu? + Bh) + 


an 
hu? + uh?) =0, 223 

re ) (223) 
in which the two equations represent conservation of mass and energy respectively. The corresponding 
Rankine-Hugoniot relations, 
fuajt _ [Zhu + un? ]* 

= ya lpayt’ 
AE ~ [hue + dna] 


give shock conditions that are quite different from (221). We have a choice of conserving either mass 
and momentum or mass and energy. In fact, two different kinds of bores are observed in practice: 
turbulent bores that conserve momentum but lose energy, and undular bores that conserve energy but 
not momentum. 


(224) 


v= 


Causality 


Once we have chosen a particular conservation form and, thus, a particular weak formulation, 
there is still the possibility of more than one solution existing if we allow shocks. As for scalar 
PDEs, there are some shock solutions that, although they satisfy the Rankine-Hugoniot 
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Figure 19: Schematics of the characteristics for two alternative shock solutions; the shock is 
shown as a heavy dashed line. 


conditions, are unphysical and should be eliminated. There are several methods for doing 
this, of which we concentrate on causality: making sure that information propagates into a 
shock, rather than out of it. 

An n-dimensional hyperbolic system has n families of characteristics, so a shock intersects 
2n of them: n from either side. If there are k outgoing characteristics, then there are (2n — k) 
characteristics going in. We also have the n Rankine-Hugoniot relations, giving a total of 
(3n — k) pieces of information on the shock. The unknowns are the n components of u, on 
either side of the shock, and the shock slope dy/d, giving a total of (2n+1). For the number 
of equations to equal the number of unknowns, we require (3n — k) = (2n + 1), that is 


k=(n-1). (225) 


This is the condition for a shock to be causal: there must be (n — 1) characteristics leaving 
the shock (and, therefore, (n + 1) going in). 

For scalar equations, n = 1 so there should be no characteristics leaving the shock, 2 
going in, as imposed previously. For two-dimensional systems, n = 2 so we need one family 
of characteristics going out of a shock and three going in. Schematics of the characteristics 
for two alternative shock solutions are shown in Figure 19. In diagram (a), three families of 
characteristics enter the shock and one leaves: this solution is causal. In diagram (b), only 
one family of characteristics propagates into the shock, so this solution is non-causal and 
should be discarded. 


Example 35 Consider a shock solution of the shallow-water equations satisfying the momentum- 
conserving Rankine-Hugoniot conditions (221). Recall that the characteristic velocities are given by 
da/dt = u+Vh. Assume the shock is moving in the positive x-direction. Then, for the solution to 
be causal, as shown in Figure 19, there should be one set of characteristics entering the shock from 
behind and two from the front. In other words, the shock should be moving faster than one of the 
characteristic speeds behind the shock and faster than both the characteristic speeds ahead of the shock. 
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Thus we obtain the following four inequalities 
u-t</S/h_>v, u—VJh_<v, uth <v, u—VWhi <2», (226) 
where, as before, v is the speed of the shock. From these it follows that 
(u_ —v)? <h_ and (u,—v)?>h4, (227) 
and the Rankine-Hugoniot condition (221) then leads to the condition 
hy < he. (228) 


Thus the height behind a bore must be greater than that ahead; otherwise the bore is not causal and 
the discontinuity cannot be maintained. That all the inequalities in (226) follow from (228), may be 
verified from the identities 


(u, —v)? = oa os, (u_—v)? = (et a (229) 


which follow from (221). 
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6 Second order semi-linear equations 


6.1 Introduction 


Now we consider second-order scalar equations of the form 
Oru Oru Oru Ou Ou 
=z + 2b — —— —,—]}. 230 
ole.) ooe + a Wee + else =f (tam > | (230) 


This is called a semilinear equation because the coefficients a, b and c are independent of 

u and its derivatives. Note that a second-order scalar equation like (230) may always be 

transformed into a first-order system by setting p = Ou/Ox, q = Ou/Oy. Where one of the 

independent variables is clearly supposed to represent time, we will also write (230) in the 
— a oO? a Ou O 

u u u u Ou 

a(x, t)— > + 2b(2, t).—— + c(x,t)—> = xz,t,u,=——,— }. 231 

(ot) + WEDS + eat) 53 =f (atu 5] (231) 

We will concentrate on three canonical examples of particular importance, namely the 

wave equation 


es = = = 0, (232a) 
the heat equation ; 

a — = = 0, (232b) 
and Laplace’s equation ; ; 

= + = =0. (232c) 


6.2 Cauchy data 


Cauchy data for the general second-order quiasi-linear PDE (230) is to specify u and its 
normal derivative on some curve [ in the (z, y) plane, i.e. 
Ou 
t= r0(s), y= yo(s), U= uo(s), an = vo(s), (233) 
where s parametrises I. This is equivalent to specifying u and both its first derivatives on IT’, 
since 


dug  Oudzg | Oudyo 
ds Oxds dy ds 


= 8 
_ Ou (Oudyo Oudzo day \? dyo s 
tole) = On (> ds Oy ds ) { ($2) r (2 ee) 


may be solved simultaneously for 0u/Ox and Ou/Oy. We may therefore replace (233) with 


(234) 


and 


r=20(s),  y=wls, w=wls), F¢=po(s), SY =anls). (236) 
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A necessary condition for the solution u(z,y) to be uniquely defined in a neighbourhood 
of I’ is for the second derivatives of u to be determined on I’. Now, if we differentiate the 
initial data along I’ and use the chain rule, we obtain 


dpo _ dao 0?u dy O7u dqg _ dz Oru _ yo Oru 
ds ds Ox? ds OxOy’ ds ds Oxdy ds Oy?’ 


(237) 


Along with the PDE (230), this gives us a system of three equations for the three second 
partial derivatives of u, and the determinant of this system is 


a 2b Cc 
dzp dyo 2 2 
a ae fa |) oe ee | (238) 
ds ds ds ds 
dzo dye 
ds ds 


A necessary condition for the Cauchy data (236) to determine uw locally is, therefore, 


dyo : dxo dyo dxo : 


6.3 Characteristics 
Definition 


As for first-order equations, we define characteristics to be curves in the (a, y) plane on which 
Cauchy data do not determine a unique solution. If such a curve is parametrised by x = x(7), 
y = y(T) then, from (239), we have 


dy ‘ dx dy dx \? _ 
a (S*) 2b An +c (=) = 0, (240) 


so the slopes of the characteristics satisfy 


dy 
ne | 
dx 


where ad? — 2b\+c=0 (241) 


Characteristics may also be defined as curves across which the second derivatives of u may 
be discontinuous, with lower derivatives continuous. It is left as an exercise to show that this 
alternative definition leads to the same equation (241) for the characteristic slopes. 


Classification 


The PDE (230) is classified according to the number of distinct real characteristics it possesses, 
and this is determined by the sign of the discriminant b? — ac. 


1. b? > ac implies that there are two distinct real characteristics and the PDE is hyperbolic. 


2. b? < ac implies that there are two complex conjugate characteristics and the PDE is 
elliptic. 


3. b? = ac implies that there is one repeated real characteristics and the PDE is parabolic. 
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Since a, 6 and ¢ are functions of x and y, the PDE may change type (i.e. hyperbolic, 
parabolic or elliptic) from one region to another. However, the type of a PDE at any given 
position is invariant under coordinate transformations. 


Example 36 Suppose we transform to new coordinates (X,Y), which are given as functions of x 
and y. Then the PDE (230) is transformed to 


Ou Oru O7u 
Aaya 2Bavay + Can = F, (242) 
where 
Ox\* Ox oax ax? 
OX OY OX OY OX OY OX OY 
Baas Ox (= oy Da mm) em Oy’ eee) 
dY\* | ay oY ay \* 
It follows that 
aX OY  axoyY\? 
2 — (72 
B? — AC = (0? — ac) ( ae oy Oy | ; (244) 


so that the sign of the discriminant is invariant. 


6.4 Canonical forms 
Hyperbolic equations 


For hyperbolic PDEs, the quadratic equation (241) has two real distinct roots, say A1(x, y) 
and A2(z,y). Now suppose the differential equations dy/dxz = ),; have first integrals given by 


dy 


ri Ai (x,y) => -€y) const, (245a) 
dy _ _ 
a A2(x, y) => (x,y) = const. (245b) 


Now we change variables from (x,y) to (€,7), using the chain rules 


Ou 0€ Ou On Ou 


be On Be On On a8) 


and so forth. The PDE (230) is clearly transformed to an equation with an analogous form, 


that is, : 
O7u O7u O-u Ou Ou 


We know that the characteristics are given by € = const and 7 = const, so the roots of the 


quadratic form 
dn \? dn dé dé \? 
2 — = 248 

a($) Pardr dr ° (248) 


(247) 
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must be dé/dr = 0 and dn/dr = 0. It follows that a = y = 0, so that (247) takes the form 


Ou Ou Ou 
= 9414, Der BH de 249 
Seay = 8 (E55) a 
This is the so-called canonical form for second-order hyperbolic PDEs. 
Example 37 For the PDE 
Ou dru 
eee 2 
Fat 7 ppp = fle (250) 
the characteristics are given by 
dy . 
(s¢) —-1=0 => yxtx= const. (251) 
So we set€=x-—y, n=x+y and, by changing variables, obtain 
Ou 1,f€E+n n-€ 
= — Pe 2 2 
aéan it ( eo ) o(§n), say (252) 
This may now be integrated directly to give the so-called D’Alembert solution 
u= ff 66.1) agan-+ ha(6) + hale). (253) 


Elliptic equations 


For elliptic equations, the roots of the characteristic equation (241) are complex conjugates, 
say 


d : 

Fp = Any) + idr(2,9). (254) 
Suppose that the integrals of these ODEs can be written in the form 

£(z,y) + in(s, y) = const, (255) 


for some functions € and 7. Then we use € and 7 as new variables. Again, the transformed 
equation must be of the form (247) but, this time, the roots of 


dn\? _, ,dndé dé\? 
- (52) ot ae dr a (¢) ad i226) 


are d€/dr + idyn/drt = 0, which implies that 8 = 0, y = a. The canonical form for elliptic 
PDEs is, therefore, 


Ou Ou Ou Ou 
ee ye a aaa 257 
Example 38 For the PDE 
Oru Oru Oru 
: 2 + 22? —— = 2 
Y Ox? Vaxay ° Oy? m eae) 
the characteristics satisfy 
dy =(1+ i)= => (l+i)2? —y? = const. (259) 
da y 
If we choose £ = x7, n = 2? — y”, then the equation transforms to 
2 2 
Oru | Oru _ 1 Ou E+ Ou (260) 
Og? On? 2E DE \ 2E(E—n) / On 
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Parabolic equations 
For parabolic PDEs, there is one repeated characteristic slope 


dy 0b 
ea aera 261 
dz a o] ( ) 
which we suppose has the solution (x,y) = const. Then any convenient linearly independent 
function €(x, y) may be chosen, and the PDE (230) transforms to the canonical form 


O7u Ou =) 


0e2 = Q (6.0, U, 0&’ on (262) 
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7 Semi-linear hyperbolic equations 


7.1 Non-Cauchy data 


As stated in Section 5, Cauchy data for second-order equations specifies both u and its normal 
derivative on a plane curve I. In general, however, there may be boundaries on which it is 
appropriate to prescribe just one BC or even none. 


Example 39 Consider a string of length L and density (i.e. mass per unit length) p, stretched to a 
tension T between the points x =0 anda =L. It may be shown that small transverse displacements 


u(x,t) satisfy the wave equation 
Oru O7u 
ae ~° Oa" en) 


where t is time and c? = T/p. Suppose we wish to determine u(x,t) in the region 0 < x < L, 


0<t<T. To do so, we need to specify the initial displacement u(x, 0) = uo(x) and the initial velocity 
Ou/Ot(x,0) = vo(a); these correspond to Cauchy data on the initial curve t =0,0<a< L. On the 
boundary curves x = 0 andx = L, we give just one boundary condition, namely that the displacement 
is zero: u(0,t) = u(L,t) =0. Finally, on the boundary t = T we give no boundary conditions at all. 


Example 39 illustrates that, in general, the number of BCs needed on any boundary is 
equal to the number of characteristic families travelling out of that boundary. Where two 
sets of characteristics travel out, a boundary is called space-like, and two conditions must be 
given. On a time-like boundary, with one characteristic family travelling in and one travelling 
out, just one condition must be given. Finally, no conditions may be given on a boundary 
where all the characteristics travel in. 


7.2. Weak solutions 


Weak solutions may be defined in a manner analogous to that used to define first-order 
PDEs. Notice that, unlike quasilinear equations, the semilinear equations considered here do 
not spontaneously form shocks. Consequently a discontinuous solution will only arise if dis- 
continuous initial data are specified. Furthermore, for semilinear equations, all discontinuities 
propagate along characteristics. 


7.3. Riemann’s method 


Riemann’s method is a way of solving linear hyperbolic equations that are stated in (i.e. have 
been transformed into) canonical form 


O7u Ou Ou 


Suppose Cauchy data are given on a curve I in the (x,y) plane, so that u, 0u/Ox and Ou/Oy 

are all known on I. The condition (239) for these data to determine the solution locally 

translates into the restriction that [ must be nowehere parallel to either the x- or y-axis. 
We define the adjoint of the differential operator £ in (264) by 


bu) + cv. (265) 
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Figure 20: The integration region D for Riemann’s method. 


The point of this definition is to make vL[u] — uL*|v] be the divergence of something. The 
general procedure for constructing the adjoint of a differential operator is (7) put all coefficients 
inside the derivatives; (7) switch the signs of all odd-order derivatives. By combining (264) 
and (265), it is readily shown that 


O 0 0 0 
v£L[u] — ul*[v] = ae (os a aw) + a (-u ~ bu) (266) 
Now we integrate over the region illustrated in Figure 20, namely the region bounded by 
T and the lines x = € and y = , where € and 77 are constants. The intersection between y = 1 
and [ is labelled A, the point where x = € intersects [is labelled B, and the point (£,7) is 
labelled P. By using (266) and Green’s theorem, we obtain 


If {RL[u] — ulL*[R]} dedy = f {Rr (5 + au) dy +u (= - on) ar} , (267) 


where R is called the Riemann function. The idea is to choose the properties of R to simplify 
(267) as much as possible. 

On the left-hand side of (267), we use the fact that Liu] = f and choose R to satisfy 
L*|R] = 0. On the right-hand side, we note that dy = 0 on AP and that dr = 0 on PB. The 
latter integral along PB we integrate by parts: 


A P 
/ Rfdady = | u (Ze oR) d+ f u(ar- 5) dy + [Ru] 
D P Ox B Oy 
B 
+ R acre dy + u of oR dxzp. 
A Oy Ox 


Now we suppose that 0R/Ox = bR on AP and that 0OR/Oy = aR on BP, to eliminate the 
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r 


Figure 21: Schematic showing the behaviour of the Riemann function. 


first two integrals. If we also choose R = 1 at the point P, we end up with 


wen II Rf dady + R(B)u(B) - fo {n(e% + au) dtu & - oR) ar} (268) 


This gives us the solution at any arbitrary point (€,7) in terms of f, u and its derivatives on 
I (t.e. the Cauchy data), and the Riemann function R, so if we can find R then the problem 
is solved in principle. 

To summarise, the properties we require of R are 


4 PR Oo 3) 
E aay 5g ak) a) t+ ch = 0, i ae af 
OR 
—=bR =); 
Ox aan (269) 
OR 
ay =aR c= E, 


R=1 (x,y) = (7). 


It may be shown that the problem (269) determines R(x, y;€,1) uniquely. Notice that R is 
independent of the function f and the boundary data applied to u: it depends only on the 
orginal differential operator £. Thus, it may be found “once and for all” and then applied to 
any boundary data and any right-hand side f. However, finding R explicitly is very difficult 
except in simple special cases. It is generally most useful for theoretical purposes, for example 
proving that u depends continuously on its initial data. 
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Finally, we note an alternative definition of R in terms of generalised functions. Suppose 
we form a domain S, in which the solution is to be found, between [ and a second curve ¥. 
Now we solve the problem 


L*|R] = 6(x — €)d(y—m) in S, 


OR (270) 
= = 0 on ¥, 
On 


for R backwards, starting from y, where 6 is the Dirac delta function. As illustrated in 
Figure 21, since R and its first derivatives are zero on y, F is identically zero outside the 
region D considered previously. Across the lines x = € and y = 7, discontinuities in R 
balance the singularity on the right-hand side of (270); in particular, 


R~H(é—«)H(n—y) as (x,y) > (€,7), (271) 


where H is the Heaviside function. It is readily verified that (270) implies that R satisfies the 
properties listed in (269). 


Example 40 The inhomogeneous wave equation in canonical form is 


Llu] = 5 = fle.y). (272) 


Here, the problem satisfied by R is 


: oR 
EN Saag Zak yon 
oF = 9 y=, 
Ox (273) 
OR 
aes r=, 


R=1 — (#,y) = (én). 


Notice that R is constant and, therefore, equal tol onx = € and ony =7. It follows that the Riemann 
function is simply given by R=1ina<€,y<yn (or R=H(€-—2)H(n-y)). Thus the solution is 


B 
u= fff tardy + w(B) - y By oe (274) 


which is equivalent to the D’Alembert solution found in Example 37. 
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8 Semi-linear elliptic equations 


8.1 Well-posed boundary data 


The canonical form for second-order elliptic PDEs is 


Oru = d?u r( Ou =) . (275) 


a OF 1th ts Be By 
The operator on the left-hand side is referred to as the Laplacian, for which the symbols V7u 
or Au are often used as shorthand. 

Recall the Cauchy—Kowalevski theorem, which states that a unique solution to (275) exists 
in a neighbourhood of an initial curve [ if u, du/Ox and Ou/Oy are specified analytic functions 
on T and satisfy the condition (239). However, the Cauchy problem is ill posed for elliptic 
equations. This manifests itself in several ways, including the following. 


e Small changes in the initial data may lead to arbitrarily large changes in the solution. 
This makes it impossible to compute meaningful solutions numerically. 


e Although a unique solution exists in a neighbourhood of I’, that neighbourhood may be 
arbitrarily small, in that a singularity may form arbitrarily close to I. 


It transpires that appropriate boundary data for (275) is to give just one boundary con- 
dition on u everywhere on a closed curve. We show below that this does indeed give a unique 
solution for Poisson’s equation, a special case of (275). This should be constrasted with 
the situation for hyperbolic equations where, as illustrated in Figure ??, either 0, 1 or 2 
boundary conditions are given at each point on a closed curve, depending on the number of 
inward-travelling characteristics. 


8.2 Uniqueness theorems for Poisson’s equation 


Poisson’s equation is a special case of (275), in which the right-hand side f depends only on 
x and y. Consider the so-called Dirichlet problem, namely Poisson’s equation for u(x, y) in 
some domain D, with u given on the boundary of D: 


Ou Ou 
5x2 + dy? = f(a, y) (2, y) ED, (276) 


u=9(z,y) (x,y) € OD. 


Vu= 


We will now show that, if a solution of (276) exists, then it is unique. 
Suppose there exist that two solutions u; and uz of (276). If dé = u, — ua, then ¢ satisfies 
the homogeneous Dirichlet problem 


(277) 


Now consider the Dirichlet integral 


V -(¢V¢) drdy = e Pees ae ge dady. (278) 
I, ee Ox} Oy \" Oy 
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This may be written in two ways (i) by expanding out the derivatives on the left-hand side, 
(ii) by using Green’s theorem: 


2 2 = o¢ 
II {oV°o+|Vo|?} dady = fos ds, (279) 


where n and s refer to the outward-pointing normal and arclength respectively on 0D. Now 
we simplify by using the properties (277) of ¢: 


If |V ol? drdy = 0. (280) 


Since the integrand is non-negative (and assumed to be continuous), it must be identically 
zero in D. It follows that ¢ is constant in D where, since ¢ = 0 on OD, that constant must 
be zero. Hence wu; = wg, so the solution to (276), if it exists, is unique. 
Next consider the Neumann problem, in which the normal derivative of u, rather than u 
itself, is specified on OD: 
Vu = f(z, y) (x,y) 2D, 


oe =g(z,y) (2,y) € OD. (281) 


By Green’s theorem, 


/ V-udady = ou ds => / f dady = fads, (282) 
D ap On D 


This is the solvability condition for (281); if (282) is not satisfied, then (281) has no solution. 
Now suppose that (281) is satisfied and let wu; and ug be two solutions of (281). Then 
o = U1 — Ug satisfies the homogeneous Neumann problem 


V7o=0 (#,y)€D, 

2 
ve =0 (a,y) € OD, os 
and it is straightforward to show (using the same approach as for the Dirichlet problem) that 
@ must therefore be constant. Thus, if the solvability condition (282) is satisfied, then the 
solution to (281), if it exists, is unique up to the addition of an arbitrary constant. 

This behaviour of the solutions of the Neumann problem is an instance of the Fredholm 
alternative. The homogeneous problem (283) has the nontrivial solution ¢ = const. Thus 
there is no solution to the inhomogeneous version (281) unless an orthogonality condition 
(282) is satisfied. If so, then the solution to (281) is nonunique since an arbitrary multiple of 
the homogeneous solution may be added. 


8.3. Maximum principle 


Suppose wu satisfies Poisson’s equation in a bounded domain D, 
Vu = f(z,y) in D, (284) 


where f > 0 in D. Then u attains its maximum value on OD. 


3 Actually, the theorem also holds if f depends on u and its derivatives, provided it is non-negative in D. 
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This was proven in the Part A DE’s course. I include the proof here for completeness, 
but it 7s non-examinable. 

The proof proceeds in two parts. First suppose that f is strictly positive in D. If u has 
an interior maximum at some point (29, yo) inside D, then the following conditions must be 
satisfied at (x0, yo): 


2 2 2 2 2 2 
Ou aM oe O"u\ (Oru > Oru Ge 26 (285) 
Ox Oy Ox? Oy? OxOy Ox?’ Oy? 
But, if f is strictly positive, (284) implies that it is impossible for both 0?u/0x? and 0?u/dy? 
to be < 0. Hence u cannot have an interior maximum within D, so it must attain its maximum 
value on OD. 


This completes the proof for the case where f > 0 in D. Now suppose that we only have 
f > 0in D, and consider the function 


€ 
v(v,y) = ula.y) +5 (@ +9"), (286) 
where € is a positive constant. Then 
Vu=fte>O inD (287) 


so, using the result just proved, v attains its maximum value on OD. Now, if the maximum 
value of u on OD is M and the maximum value of (x? + y?) on OD is R?, then the maximum 
value of v on OD (and thus throughout D) is M + (€/4)R?. In other words, the inequality 


uto (2+) =vSM+ 2 (288) 


holds for all (x,y) € D. Letting « > 0, we see that u < M throughout D, i.e. that u attains 
its maximum value on OD. 

It clearly follows (by using the above result with u replaced by —u) that, if f < 0 in D, 
then u attains its minimum value on OD. In the case f = 0, u therefore attains both its 
maximum and minimum values on OV. This is an important property of Laplace’s equation 
V’u = 0: in any closed region D, u is bounded between its maximum and minimum values 
on OD. It may also be used to prove uniqueness of solutions of Poisson’s equation. In the 
homogeneous Dirichlet problem (277), ¢ must take its maximum and minimum values on 0D. 
But ¢=0 on OD and must, therefore, be identically zero. 


8.4 Green’s functions 


Green’s functions play the same role for linear elliptic PDEs that Riemann functions play for 
hyperbolic PDEs. Consider the Dirichlet problem 


Vu = f(z,y) (x,y) € D, 
u= g(x,y) (x,y) € a 0282) 


Notice that the Laplacian is a self-adjoint operator, so we consider the integral 


/I {uV’G — GV7u} dady = /I V -{uVG—GVu} dady 
Ss S 


OG Ou 
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where n and s refer to the outward-pointing normal and arclength on OD, and G is the 
Green’s function. If we choose G to satisfy 


WG = (a — £)6(y — in D 
(x—€)d(y—n) inD, (291) 
G=0 on OD, 
where 6 is the Dirac delta-function, then we have 
u(&,7) = ff Gf dady + fo gee ds. (292) 


If the Green’s function G can be found from (291), then (292) gives the solution at any point 
(€,7) € D in terms of G and the given functions f and g. 

Now we have to find a function G whose Laplacian is equal to a delta function. Thus 
G should satisfy V?G = 0 when (x,y) 4 (€,7) and should have some sort of singularity as 
(x,y) > (€,7), and it turns out that the correct singularity is 


1 
aoe log|(x, y) — (€,)| + O() as (x,y) > (&,7). (293) 
It is readily verified that log] (a, y) — (&, n)| satisfies Laplace’s equation for all (x,y) 4 (€,7). 


Hence 
/ i oV’G dady = / | VG dady, (294) 
D € 


where S, is a disc of radius € centred at (€,7) and (x,y) is any suitably differentiable test 
function. Recalling that V?G is a distribution, we evaluate the right-hand side as 


Og 
VG dady = GV? ¢odaxd (05 a - ot ) as 295 
Ife xdy iv ) ray +p p a (295) 


Now introducing local polar coordinates « = €+rcosé, y=7+rsin9, letting « > 0 and 
using the asymptotic behaviour (293), we find that 


27 € 
| VG dady ~ | i am log(r)V2or drdé 
oA 0 Jo 27 


27 al an 1 Od 
+ / 5 (én) ara / = loal) Ped. (296) 


The first and third integrals on the right-hand side vanish as € > 0, so we are left with 


/I oV’G dady > $(E,7) as e300, (297) 


and hence (294) reduces to 
i, i. oV°G dady = 4(€,7). (298) 
D 


Thus V?G is indeed equal to a delta function. 
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An alternative statement of (291) is, therefore, 
VG =0 in D\(§,n), 


G~ 5 los|(2,y) — (En) as (x,y) + (€,). 


The trick to obtaining G is thus to find a function that satisfies Laplace’s equation inside D 
and is equal to —1/(27) log|(x, y) — (€,7)| on OD. 


Example 41 Green’s function for a half-plane 
For the problem 
Veu=f(t,y) yy >9, 


u= g(x) y = 0, (300) 
u—>0 yorou , 
the Green’s function satisfies 
VG =0 y>0, (x,y) 4 (En), 
at y=0, (301) 


1 
Gw x, esl (a) — (€,7) ’ 


This problem is readily solved using the method of images. Consider an image singularity at the 
point (€,—n), that is the function 1/(27) log| (a, y) —(&, —n)|. This function clearly satisfies Laplace’s 
equation away from the singularity (which is outside the half-plane in which G is to be defined). It is 
also equal to 1/(27) log| (zx, y) — (, n)| on the line y =0. Hence the Green’s function for this problem 
1s 


(x,y) > (§,1). 


Ga, y3€,n) = 5 los (0,9) —(€,)| - 5 loa|(0,) — (€,-n)| 
1 (2 —£)? +(y—1)? 
~ 78 (i ~OF+ (yt i) : oe 


Example 42 Green’s function for a circle 
For a Dirichlet problem on a circular disc of radius a, the Green’s function satisfies 


wea =0 x oe y? < a’, (x,y) # (,1); 
G=0 a? +y* =a’, (303) 
1 
TT 


The point (€,7) inside the disc has a corresponding image point (€',1') outside the disc, defined by 


()=2%()) om 


This image point has the property that, for any point (x,y) on the circle x? + y? =a’, the ratio 


(x,y) — (&)| 
|(a,y) _ (é,77')| 


is constant, and equal to [(é, n)|/a. The Green’s function is therefore given by 


(x,y) — (€,')| fe) 


> (wy) > (én). 


a 


Gla, yn) = 5 log|(x,9) (é,n)| ee ( 


2 gs ((a — €#) + (y—n)?)a? 
~ de 8 \ + P(E +) — 2aPak — 2a2yn + af |” 


(305) 
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Notice that, unlike the Riemann function for hyperbolic PDEs, the Green’s function is 
dependent on the shape of the region D. Unless D is a simple shape, it is usually difficult 
to find G explicitly. One useful approach is to use conformal mapping to transform D into a 
simple domain such as a disc — see below 

It is also possible to define a Green’s function for Neumann problems, but this is more 
difficult since the solvability condition (282) must be enforced. This may be achieved by 
defining G to satisfy 

VG =6(x4—€)6(y—n)-C 1inD, 


306 
Ss =0 on OD, a“ 
On 


where C' is chosen to satisfy the solvability condition. Integration of (306) over D leads to 
the condition C = 1/A, where A is the area of D. Then, by an argument analogous to that 
used for the Dirichlet problem, we obtain the solution as 


u(é,n) =a+/[ Grardy- gas. (307) 


1 
= a [f waeas. (308) 


which remains indeterminate (recall that the solution of the Neumann problem is only deter- 
mined up to the addition of an arbitrary constant). 


Here, w is the average value of u, 


md] 


8.5 Complex variable methods 


There is a close connection between Laplace’s equation and functions of a complex variable 
z =x + iy. If we transform variables from (x,y) to (z,Z), where 


z=ax+iy, Z=a2-iy, (309) 


then we find 
Oru Ou _ Ou 


= =0 310 
Ox? uy Oy? 0z0z , ty) 

and it follows that the general solution is 
u = hy(z) + ho(Z), (311) 


for some arbitrary functions h, and hg. If we require u to be a real-valued function of « and 
y, then we must have hg = hy, so the general real-valued solution of Laplace’s equation is 


i Rf (z)| (312) 


(where f = 2h;). Solutions of Laplace’s equation may sometimes by found by spotting a 
function f(z) whose real part is equal to a given function on a given curve. 


Example 43 Find a solution u(x, y) of Laplace’s equation that is equal to |x| on y = 0. 
The complex-valued function 


f(z) =2+ 7 clog2 (313) 
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-3 


Figure 22: Contour plots of the solution u given by (318) versus € = R(¢) and 7 = S(C) (left 
plot); x and y (right plot). The contour values are u = —0.9, —0.8,--- ,0.8, 0.9. 


may be split into its real and imaginary parts as follows: 


fe)={(1 = \e Hiogr} +i{ (1 | y+ wer} (314) 


where (r,@) are plane polar coordinates, i.e. x = rcosé, y=rsin@. Notice that, on y = 0, R(f) ts 
equal to x when 0 = 0 and —x when 0 = 7; in other words, R(f) is equal to |x| when y =0. A suitable 
solution is, therefore, 


u(x, y) = R[ f(z)] = (1 _ ? tan™*(u/2)) Lu log (a? + y’) ; (315) 

Another useful technique from complex variable theory is conformal mapping. This is a 

mapping from one complex variable z = x + iy to another ¢ = € + in, given by a functional 

relation ¢ = g(z), where g is analytic and g'(z) #4 0. Conformal mapping may be used to 

transform a complicated solution domain to a simple one (such as a half-plane or a disc). 

This technique works because Laplace’s equation is invariant under conformal mapping; see 
a textbook on complex analysis e.g. Priestley’ for more details. 


Example 44 Consider the solution of Laplace’s equation for u(x,y) with u > 0 as 27 + y? > 00 
and u = 2x on the line segment y= 0, -l1<a<1. This line segment is the image of the unit circle 
\¢| =1 under the conformal map 


btiy=2=5(¢+3) > C=z4+vV22-1. (316) 


The boundary conditions are mapped to 


u=R(C) on |¢| = 1, ur0 asG>o, (317) 


4H. A. PRIESTLEY, 1990 Introduction to Complex Analysis, revised edition. Oxford University Press. 
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Figure 23: Schematic of the conformal transformation w = g(z) mapping D to the unit disc 
and ¢ to the origin. 


and the solution may then easily be spotted. On the unit circle, we have ¢ = ¢~! and, hence, 
R(C-1) = R(C). A suitable solution for u is thus 
w=R(C-1) = (2- V2?=1). (318) 


Contour plots of u in the ¢- and z-planes are shown in Figure 22. 


Conformal mapping may be used to find the Green’s function for any domain that can 
be conformally mapped onto the unit disc. If we write G = oral f (z)], then the problem (299) 
implies the following properties for /: 


K(f) =0 on OD, 
1 (319) 
f ~ =< log(z—¢), as 2 => ¢, 
27 


where ¢ = € + in. Now suppose we write w = g(z) where the conformal map g transforms 
D to the unit disc and the point ¢ to the origin, as illustrated in Figure 23. In the w-plane, 
(319) becomes 
K(f) =0 on |w| = 1, 
1 (320) 
f ~ — logy, as w > 0, 
27 


and the solution is simply f = (1/27) logw. Hence, if we can find the appropriate conformal 
mapping g, then the Green’s function is 


Gx = log|9(2)]. (321) 


Example 45 Green’s function for a circle (again) 
The Mobius transformation 
w= 9(z) = Se-9) (322) 


~ Cz—a2 
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Figure 24: Schematic of the conformal transformation w = g(z) mapping the strip x > 0, 
0<y<a to the unit disc and ¢ to the origin. 


maps the disc |z| < a onto the unit disc |w| < 1 and maps the point z = ¢ tow = 0. The Green’s 
function for a disc of radius a is, therefore, 


a(z — 6) 


Cz —a? 


G = — log 


2 
= , (323) 


and it is straightforward to show that this reproduces (305). 


Example 46 Green’s function for a strip 
The semi-infinite strip x > 0,0 < y <a is mapped to the unit disc (as illustrated in Figure 24) by the 
transformation 

cosh(7z/a) — cosh(m¢/a) 
cosh(mz/a) — cosh(1¢/a)’ 


w=g(z)= (324) 


which also maps z =¢ tow =0. The Green’s function for this shape is then G = 1/(27) log |g]. 
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9 Semi-linear parabolic equations 


9.1 General properties 


In parabolic PDEs one independent variable usually represents time, and we use x and t 
as independent variables instead of x and y. The canonical form for second-order parabolic 
equations is 


Oru Ou Ou 
ee | eee 325 
Ox? («. Oe 7) N22 
and specific examples include the heat equation, or diffusion equation 
du Ou 
== 326 
Ot Ox?’ ee) 
the reaction-diffusion equation 
Ou OPu 
— = —, t 326b 
Ot Ox2 + f(x, ih) ( ) 
and the reaction-convection-diffusion equation 
Ou du Ou 


Another example is the Black-Scholes equation in mathematical finance. All these equations 
have the repeated characteristic t = const. They also have the general property that any initial 
singularities in u at t = 0 are instantaneously smoothed out when ¢ is positive. Conversely, 
backward diffusion equations such as 
2 
ay = i, (326d) 
Ot Ox? 


are ill posed, with smooth initial data giving rise to finite-time singularities in general. 


9.2 Well posed boundary data 


Typical boundary data for a diffusion equation are: 


e An initial condition for u at t = 0 


e One boundary condition on each of two curves Cy and C2 in the (z,t)-plane that are 
nowhere parallel to the x-axis. 


Example 47 The heat equation 
du Ou 
—=s5 (327) 
Ot = Ox? 
is a simple model for the temperature u(x,t) in a uniform bar of conductive material, where x is 
position and t is time. Suppose the bar is of length L, its initial temperature is ug(x) and its ends are 


kept at zero temperature. Then the boundary and initial conditions are 


U = Uo(Z) t=0, (328a) 
w=0 #=10, (328b) 
u=0 c=L. (328c) 


If, instead of being held at constant temperature, an end is insulated, then the Dirichlet boundary 
condition u = 0 is replaced by the Neumann boundary condition Ou/Ox = 0. Also, the boundary 
conditions at x = 0 and x = L may in general be given on moving boundaries, say x = x(t) and 
L = 2(t). 
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9.3. Uniqueness theorem 


We will now show that the problem 


2 
oe = SE 4 Flat) DSU ty) ae eo 
ti =up(2) t=0,-94(t) < @ < walt), (329) 
= @i{t) t> 0) 2 =e), 
U = go(t) ¢>0, 2 = olf), 


has at most one solution. Suppose there are two solutions u, and ug; then ¢ = uj; — ug satisfies 
the homogeneous problem 


2 
oo 38 nim eee ety 
b= t=0,2i() <2 <a, (330) 
¢=0 LS 0, esate), 
¢=0 tS 0 e =a). 


Now, Leibnitz’ rule gives 


ie z2() | do daz 5 dz 15 
fe = ogo? dag _ dey 
dt [, o dx i: °>: dx + oF o* (xo, t) re b°(x1,t) 


x(t) Od dg x2 x2(t) ad 2 
= 2/9 eas 
i Paget | 05 i : (=) 2 


tL 
xr2(t) 2 
mee i; 2 (3) ie (331) 
x(t) Ox 
x(t) 
/ ob? dx 
x1(t) 


is clearly non-negative, is zero when t = 0 and is a non-increasing function of t. From this we 
can conclude that it must be identically zero and, hence, that ¢ = 0. 


So the integral 


9.4 Maximum value theorem 


Suppose wu satisfies the inhomogeneous heat equation 


Ou Oru 
Ot Ox? 
in a region D bounded by the lines t = 0, t = 7 and two non-intersecting smooth curves C 
and C2 that are nowhere parallel to the x-axis. Suppose also that f <0 in D. Then u takes 
its maximum value either on t = 0 or on one of the curves C}; or C4. 
This was proven in Part A DEs. The proof included here is non-examinable. 
The proof is similar to that for Poisson’s equation. We start by supposing that f is strictly 
negative in D. At an internal maximum inside D, u must satisfy 
Ou Ou | Oru Oru 


cee —~ < =z <0. 
Oz Ot m Ox? <0, Ot? a 8) 


+ f(a, t) (332) 


B5.2 Applied Partial Differential Equations 70 


On the other hand, if u has a maximum at a point on t = 7, then it must satisfy 


Ou Ou Ou 

—=0 — >0 =z <0 334 

Ox , a Ox? — i384) 
there. With f assumed to be negative, both of these lead to contradictions, and it follows 
that u must take its maximum value somewhere on OD but not on t = T. 


Now, if the inequality is not strict, i.e. f <0, then define 
v(a,t) = u(x,t) + 5a, (335) 


where € is a positive constant. Then v satisfies 


az t fet) —€ <0 (336) 


so that v takes its maximum value on either t = 0, Cy or Co. Thus, if the maximum value of 
u over these three portions of 0D is M, and the maximum value of |x| on C, and C2 is L, 
then 2 

usv< su + M. (337) 


Now we let « + 0 and conclude that u < M, i.e. u takes its maximum value on C) U C2 U {t = 0}. 

If f > 0 in D, then a similar argument shows that u attains its minimum value on either 
C1, Co or t = 0. Thus, for the homogeneous equation in which f = 0, u attains both its 
maximum and its minimum values on Cy U C2 U {t = 0}. This property may be used in a 
uniqueness proof just as for Poisson’s equation. 


9.5 Green’s functions 


Consider the problem 


2 
lu) = SY - 9 = F(a, (S071) ae =< ah), 
u=wl2)  t=0, n(t)<a<aalt), (338) 
u= gilt) t>0, r=2\(t), 
u = g(t) t>0,e= a(t): 


If the adjoint differential operator is defined by 


. aG eG 
LNG] =—B, — Bar (339) 


then Green’s theorem leads to 


[f $eetul - wertay} ava = [fff Feu | o (a ox) h dadt 
=$ 4($ oo) dt Gude} (340) 
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Increasing t - ¢ 


Figure 25: The Green’s function (345) plotted versus x when € = QO, for 
(7 — t) = 0.2, 0.4, 0.6, 0.8, 1.0. 


where D is the region bounded by t = 0, t= 7, « = x1(t) and x = x(t). So, if we choose the 
Green’s function G to satisfy 


; 0G 0G 
—L*|G| = At at = 0 in D, 
G=0 x= x(t), (341) 
G=0 c= x(t), 


G = 6(x - €) =F 
where 6 is the Dirac delta-function, then equation (340) reduces to 


x2(0) 


oe / I Ce flaw dedi / Ceoe Auta 


1(0) 


+ [ nO Fed.tsrat— [oto Zteat.tg rat. (342) 


If G can be found from (341), then this gives an explicit formula for the solution u at 
any point (€,7) € D. Notice that G satisfies a backward diffusion equation and thus is solved 
backwards in time, starting at t = 7. As for elliptic equations, actually finding the Green’s 
function is difficult except in a few simple cases. 


Example 48 The simplest case is when the heat equation is to be solved on the whole real line, i.e. 


2 
ou F8 = F(a,1) t>0, —c<£< 00, 


U = uo(x) t=0, -w<4<ow, (348) 
u—>0 t >0, x — too. 
Then the Green’s function satisfies 
0G @G ~ . 
—{.- =~ =0 t>0, -w<%<oM, 
Ot Ox? 
(344) 


0 
G = 6(2) t=0, —co<%< 00, 
G0 £>0 


, &£— +00, 
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Increasing t - ¢ 


Figure 26: The Green’s function (349) plotted versus x when € = 1, for 
(7 —t) =0.2, 0.4, 0.6, 0.8, 1.0. 


where & = «—€ andt=1—t. This problem corresponds to the heat flow due to an initial “hot spot” 
concentrated at x = 0. The solution may be found, for example, by taking a Fourier transform in x 


(see Example ??): 
(x — £)? ) . (345) 


os 0 (— A(7 — t) 


whose behaviour is illustrated in Figure 25. Thus the solution of (343) is 


usr) =f Mew (-S (@— 2) a +f Dawes »( sr) dxdt. (346) 


Example 49 Nezt consider the heat equation on a half-line: 


oe FY = f(a, t>0,0<4%<0o, 
U = Uo(x) t=0,0<4<a, (347) 
u = g(t) t>0, «=0, 
u—>0 t>0, r7 ~w, 


ee t<T,0<2%<o 

Ot Ox? , é 
G = (a — ) t=7T,0<24<~a, (348) 
G= t<7T, x=0, 
G0 t<7T, £00. 


Now, the Green’s function (345) found previously satisfies the heat equation and has the right 
singular behaviour as (x,t) > (7), but it is not zero when x = 0. However, we can enforce this 
boundary condition by putting an image singularity at the point (—€,T), that is 


~ 9 — exp ( a 2 => exp ( wi, ee 


which is illustrated in Figure 26. This Green’s function can be used to reproduce the solution (?7), 
(??) of Example ??. 
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9.6 Similarity solutions 


If the PDE and boundary conditions have a certain symmetry, then a similarity solution may 
be sought in the form 


u(x, y) =a x f(n), y= o (350) 


and the PDE for u becomes an ODE for f. 
Consider the heat equation 


Ou Ou 


and introduce the transformation 
€=e2, tr=et and w(f,7) = ule ~,e°r), € ER. 


This change of variables gives 


Ou _.0wOr _ 4 .dw Ou_ _,dw0€ _ ac OW 


O° Orat © Or’ On ° OB Ox OE’ 


Oru 
Or 
The heat equation transforms to give 


2 2 
boc Ow _ ga-cFw or b—c Ow c2a-v9 Ww 
OT 0&2 OT O€? 


a-cO'w OF _ 24-0" 
02 Ox = OEE 


and =e 


and is invariant under the transformation Ve if b = 2a. Thus, if u solves the equation at x,t 
then w = e~°u solve the equation at « = e~%€, t = €-°r. 
We can construct groupings of independent variables which are invariant under this trans- 


formation: 
E en x 


a/b _ (ebt)a/6 = $a/b 


which defines the dimensionless similarity variable (x,t) = x/\/t since b = 2a. [Note: 
n> co ifx—>coort>0and7n=0 if x=0.] Also, 


W . sete |e () 
zed — (ebt)e/b — ye/b ~ ONE 


suggests that we seek a solution of the heat equation of the form u(a,t) = t°/?¢v(n). Indeed, 
since the heat equation is invariant under the transformation, we anticipate that the solution 
will also be invariant under the transformation. The partial derivatives transform as follows: 


Ou C ie/2a— c/2a @) 1 c/2a— c 
em 5gt tUn) + PAu! (nt = St (Ev(n) — nv"), 


since on = —a/(2t3/2) = —n/2t, and 


Ou c/2a 0 c/2a— Pu c/2a— 
an = 4e/2 v(m) — 4¢/2 1/2! a =t /2 Ly(n). 


n); Ox 
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Then the heat equation reduces to an ODE 


#1 (2u"(n) + nu!'(n) — yv(n)) = 0, (352) 


with y = c/a such that u(a,t) = t7/2u(n) and n = x/\/t. We conclude that we may be able 
to solve the heat equation using (352) if we can write the boundary and initial conditions for 
u(x,t) as conditions on v and 7. 


Example 50 Consider the heat equation 
ou _ au 
Ot Ox?’ 


subject to u(x,0) = u(co,t) = 0, u(0,t) =t™ for some constant m. 
We seek a solution of the form 


x 


u=t™ ; =—, 353 
fim), 7 Vi (353) 
for some function f. If this form of u is substituted into the heat equation, f is found to satisfy 
d?f ndf 
+" _=mf= 4 
dyn? = 2dn Foo (02) 


with the boundary conditions f(0) = 1, f > 0 as 7 > co. The simplest case is m = 0, when the 
solution is . so 


= — e* ds = erfe 7/2), 355 
Taw (n/2) (355) 


f(n) 
the complementary error function. 
Example 51 Consider the problem 

Ut =Uge, ond<ax<ow,t>QdJ, 
with u(z,0)=0 for0<a<m, 
6) 
aay atx =0,t > 0, 
Ox 
andu>0O asx—>oo, Vt>0. 
As above, we seek solutions of the form 


u(z,t)=t"f(m), n= 7 


for some function f. By substituting this form of u in the heat equation we deduce that f satisfies 


ay nds 


i 3a (356) 
Now 5 : 
OU _ ym-1/2 gr ou _ ym—1/2 pry 
see en) SEP =e NP 0) = = 


x=0 
Since q does not depend ont, we require m = 1/2 so that u(x,t) = Vtf(n) where f solves 
aie ie 


a ee ees = , / = ~— 
ap 5 ° 0, with f’(0) q and f +0.as7 > o. 
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75 
Since f = 1 solves the ODE, we seek a solution of the form f(n) = ng(n) where 
g’ _ _(2+n?/2)_ 2 
g! n n 2 
Integrate with respect to 7 to get 
kxe 1/4 7 e@/4 
g= re > a(n) = ky +h | 92 


eo /4 
=> 9(n) =Kot ki r + 


” 2 2 a 2 
| = a) and f(n) = kon + ki ( es +nf e (a0) ‘ 
0 0 


The constants of integration Kg and Kk, are determined from the initial conditions. 


f=nfi+g => f'(0)=—-¢=9(0) =o. 
Also, 


fin)~n (1 + nf *1*9) =0 > = -Kov7, 
0 


since ia e9/4d9 = 24. e-" dr =n. The solution of the equation becomes (eventually) 


f(n) =-- 


t 2 ue re 2 
u(x,t = 20/4 e” wit f e’ dr |. 
( ) T ( vi x /2/t 
Example 52 For the telegraph equation 


Ou 
(l= 5g tea flo) (357) 
the Riemann function satisfies 
OP7R 
*lu] = ——+R=0 358 
f= aay <i y<n, (358) 
with, as in Example 40, R=1 onx=€ and ony=n. Alternatively, ift=€—-—2z, y=n-y, then 
O?R 
R=0 
[i 


R= 


R= y =0. 


Notice that this problem is invariant if the two independent variables x and y are swapped. It follows 
that a similarity solution of (359) exists in the form 


R= F(s), 
and F satisfies 


where s =2,/ry, 


(360) 
dF ildF 
—+-—4+F=0 0 
ds? s ds , sae 
with F(0) = 1 
telegraph equation is 


(361) 
The solution is the Bessel function F = Jo(s),° so the Riemann function for the 


na {ne (E=a)(n-0)) @<& y<m, 
0 


otherwise. 


(362) 


the other linearly independent solution to this second-order ODE is Yo(s), which is singular as s > 0 


